Circulant SSA

x <- 1:100
w <- 1/20
sin(2*pi*w * x + 1/2) + cos(2*pi*w* x + 10*3.3412) |> plot()

numeric(0)

Вспомогательные функции

library(Rssa)

Загрузка требуемого пакета: svd
Загрузка требуемого пакета: forecast

Присоединяю пакет: ‘Rssa’

Следующий объект скрыт от ‘package:stats’:

    decompose

library(signal)

Присоединяю пакет: ‘signal’

Следующий объект скрыт от ‘package:Rssa’:

    roots

Следующие объекты скрыты от ‘package:stats’:

    filter, poly

library(gsignal)

Registered S3 methods overwritten by 'gsignal':
  method         from  
  plot.grpdelay  signal
  plot.specgram  signal
  print.freqs    signal
  print.freqz    signal
  print.grpdelay signal
  print.impz     signal
  print.specgram signal

Присоединяю пакет: ‘gsignal’

Следующий объект скрыт от ‘package:signal’:

    Arma, as.Arma, as.Zpg, bartlett, bilinear, blackman, boxcar, butter, buttord, cheb1ord, chebwin, cheby1, cheby2, chirp, conv,
    decimate, ellip, ellipord, fftfilt, filter, filtfilt, fir1, fir2, flattopwin, freqs, freqs_plot, freqz, freqz_plot, gausswin,
    grpdelay, hamming, hanning, ifft, impz, interp, kaiser, kaiserord, levinson, Ma, medfilt1, poly, remez, resample, sftrans, sgolay,
    sgolayfilt, specgram, triang, unwrap, Zpg, zplane

Следующие объекты скрыты от ‘package:stats’:

    filter, gaussian, poly

source("eossa_new.R")


dftmtx <- function(n) {
  y <- stats::mvfft(diag(1, n))
  y
}

diag_averaging <- function(A){
  B <- A[nrow(A):1, ] |> Re()
  lapply(split(B, -(row(B) - col(B)) ), mean) |> as.numeric()
}

shift_vector <- function(vec) {
  last_element <- tail(vec, 1)
  vec <- vec[-length(vec)]
  shifted_vec <- c(last_element, vec)
  return(shifted_vec)
}

extend <- function(x, H){
  # Вычисление коэффициентов AR модели для дифференцированного ряда
  N <- length(x)
  p <- floor(N / 3)
  dx <- diff(x)
  # A <- ar(dx, order.max = p, method = "yule-walker")$ar
  A <- aryule(dx, p)$a
  
  # Правое расширение
  y <- x
  dy <- diff(y)
  er <- signal::filter(A, 1, dy)
  dy <- signal::filter(1, A, c(er, rep(0, H)))
  y <- y[1] + c(0, cumsum(dy))
  
  # Левое расширение
  y <- rev(y)
  dy <- diff(y)
  er <- signal::filter(A,1,dy)
  dy <- signal::filter(1,A,c(er, rep(0, H)))
  y <- y[1] + c(0, cumsum(dy))
  
  # Расширенный ряд
  xe <- rev(y)
  
  # Вывод результатов
  xe 
}

CiSSA

Подаётся на вход временной ряд, длина окна (если её нет, то она равна длине ряда + 1 пополам) и информация о том, нужно ли расширить ряд. Расширять ряд стоит при стохастическом тренде (Autoregressive extension (default). It is indicated for stationary and stochastic trend time series as well). Реализовано только Autoregressive extension.

На выходе список выдаётся список list(t_series, importance).
t_series — матрица, по столбцам которой располагаются временные ряды, отвечающие за частоты (i-1)/L, где i — номер столбца, L — длина окна.
importance — вектор, отвечающий за значимость i-ого временного ряда в разлолжении. Чем больше значение, тем больший вклад внёс i-тый временной ряд.

circulant_SSA <- function(ts, L = NULL, extend_flag = FALSE){
  time_series <- ts
  # Construct trajectory matrix
  N <- length(time_series)
  if (is.null(L)){
    L <- (N + 1)%/%2
  }
  # Проверка на расширения ряда
  if (extend_flag == FALSE){
    H <- 0
    time_series <- ts
  }
  else{
    H <- L
    time_series <- extend(ts, H)
  }
  
  X <- hankel(time_series, L)
  
  # Number of symmetric frequency pairs around 1/2
  if (L %% 2) {
    nf2 <- (L + 1) / 2 - 1
  } else {
    nf2 <- L / 2 - 1
  }
  
  # Number of frequencies <= 1/2
  nft <- nf2 + abs((L %% 2) - 2)
  
  # Decomposition
  # Estimate autocovariance     OK
  autocov <- numeric(L)
  for (m in 0:(L-1)){
    autocov[[m+1]] <- sum(time_series[1:(N-m)] * time_series[(1+m):N]) / (N-m)
  }
  
  # First row of circulant matrix
  circ_first_row <- numeric(L)
  for (m in 0:(L-1)){
    circ_first_row[[m+1]] <- (L-m)/L * autocov[[m+1]] + (m)/L * autocov[[L-m]]
  }
  
  # Build circulant matrix
  S_C <- matrix(circ_first_row, nrow = 1)
  shifted_vector <- circ_first_row
  for (i in 2:(L)) {
    shifted_vector <- shift_vector(shifted_vector)
    # S_C <- rbind(S_C, as.vector(shifted_vector))
    S_C <- rbind(as.vector(shifted_vector), S_C)
  }
  
  # Eigenvectors of circulant matrix (unitary base)
  U <- dftmtx(L)/sqrt(L)
  
  # Real eigenvectors (orthonormal base)
  U[, 1] <- Re(U[, 1])
  for (k in 1:nf2) {
    u_k <- U[, k + 1]
    U[, k + 1] <- sqrt(2) * Re(u_k)
    U[, L + 2 - (k + 1)] <- sqrt(2) * Im(u_k)
  }
  if (L %% 2 != 0) {
    U[, nft] <- Re(U[, nft])
  }
  
  # Eigenvalues of circulant matrix: estimated power spectral density
  psd <- abs(diag(t(U) %*% S_C %*% U))
  
  # Principal components
  W <- t(U) %*% X
  # Reconstruction
  # Elementary reconstructed series
  R <- matrix(0, nrow = N+2*H, ncol = L)
  for (k in 1:L) {
    R[, k] <- U[ ,k] %*% t(W[k, ]) |> diag_averaging()
  }
  
  # Grouping by frequency
  # Elementary reconstructed series by frequency
  Z <- matrix(0, nrow = N+2*H, ncol = nft)
  Z[, 1] <- R[, 1]
  # Importance of component
  imp <- numeric(nft)
  lambda_sm <- sum(psd)
  imp[1] <- psd[1]/lambda_sm
  for (k in 1:nf2) {
    Z[, k + 1] <- R[, k + 1] + R[, L + 2 - (k + 1)]
    imp[k+1] <- (psd[k+1] + psd[ L + 2 - (k + 1)])/lambda_sm
  }
  if (L %% 2 != 0) {
    Z[, nft] <- R[, nft]
    imp[nft] <- psd[nft] / lambda_sm
  }
  
  list(t_series = Z[(H+1):(N+H),],
       importance = imp,
       freq = (0:(length(imp) -1))/L
       )
}

# groups - list of frequencies
grouping_cissa <- function(cissa_res, groups){
  freq <- cissa_res$freq
  t_series <- cissa_res$t_series
  
  
  residuals <- 0
  result <- setNames(as.list(rep(0, length(groups))), names(groups))
  result_freqs <- list()
  for (i in 1:length(cissa_res$freq)){
    flag <- FALSE
    for (name in names(groups)){
      if (groups[[name]][1] <= freq[i] & freq[i] <= groups[[name]][2]){
        flag <- TRUE
        result[[name]] <- result[[name]] + t_series[, i]
        result_freqs[[name]] <- c(result_freqs[[name]], freq[i])
      }
    }
    
    if (flag == FALSE){
      residuals <- residuals + t_series[, i]
    }
  }
  
  result[["residuals"]] <- residuals
  
  return(
    list(
      t_series = result,
      freqs_by_group = result_freqs
    )
  )
  result
}

generate_ts <- function(func, n=1e3, ...){
  1:n |> func(...) |> ts()
}

f_cos <- function(x, A = 1, omega = 1/4, phi = 0){
  f_exp_mod_harm_series(x, A, alpha = 0, omega = omega, phi = phi)
}

f_sin <- function(x, A = 1, omega = 1/4, phi = 3*pi/2){
  f_exp_mod_harm_series(x, A, alpha = 0, omega = omega, phi = phi)
}

f_exp <- function(x, A = 1, alpha = 1){
  A * exp(alpha * x)
}

f_exp_cos <- function(x, A = 1, alpha = 1, omega = 1/4, phi = 0){
  f_exp_mod_harm_series(x, A, alpha, omega, phi)
}

f_const <- function(x, C = 0){
  rep(C, length(x))
}

f_exp_mod_harm_series <- function(x, A = 1, alpha = 1, omega = 1/4, phi = 0){
  A*exp(alpha*x)*cos(2*pi*omega*x + phi)
}

f_linear <- function(x, a = 1, b = 0){
  a*x + b
}
mse <- function(f_true, f_reconstructed){
   mean((f_true - f_reconstructed)^2) 
}

Ошибка при Lw in N, Kw not in N

n <- 96*2+5
L <- 96
eps <- 1/97
f_sum <- function(x){
  f_const(x, C = 1) + f_cos(x, omega = 1/12) 
}


f_const |> generate_ts(n, C = 1) |>
  plot(col = "green", ylim = c(-1, 2), ylab = "f_n")
f_cos |>
  generate_ts(n, omega = 1/12) |>
  lines(col="green")

f_sum |> generate_ts(n) |> lines(lwd = 3, col='red')
f_n <- f_sum(1:n)



c <- circulant_SSA(f_n, L = 96, extend_flag = FALSE)
r <- grouping_cissa(c,
               groups = list(
                 trend = c(0, eps),
                 sesonal = c(1/12-eps, 1/12 + eps)
               )
               )$t_series

f_C <- f_const |> generate_ts(n, C = 1)
f_c <- f_cos |> generate_ts(n, omega = 1/12)
print("Ошибки при CiSSA")

[1] "Ошибки при CiSSA"

print(paste("Ошибка при вычислении C = 1: ", mse(f_C, r$trend) |> format(scientific = TRUE, digits = 2) ))

[1] "Ошибка при вычислении C = 1:  3.2e-31"

print(paste("Ошибка при вычислении cos(pi/12): ", mse(f_c, r$sesonal) |> format(scientific = TRUE, digits = 2) ))

[1] "Ошибка при вычислении cos(pi/12):  5.1e-30"

lines(1:n, r$trend, col="blue")

lines(1:n, r$sesonal, col="blue")

f_const |> generate_ts(n, C = 1) |>
  plot(col = "green", ylim = c(-1, 2), ylab = "f_n")
f_cos |>
  generate_ts(n, omega = 1/12) |>
  lines(col="green")

f_sum |> generate_ts(n) |> lines(lwd = 3, col='red')
f_n <- f_sum(1:n)

s <- ssa(f_n, L = 96)
r <- reconstruct(s, groups=list(
  trend = 1,
  sesonal = 2:3
))


print("Ошибки при SSA")

[1] "Ошибки при SSA"

print(paste("Ошибка при вычислении C = 1: ", mse(f_C, r$trend) |> format(scientific = TRUE, digits = 2)  ))

[1] "Ошибка при вычислении C = 1:  9.6e-05"

print(paste("Ошибка при вычислении cos(pi/12): ", mse(f_c, r$sesonal) |> format(scientific = TRUE, digits = 2)))

[1] "Ошибка при вычислении cos(pi/12):  9.6e-05"

lines(1:n, r$trend)

lines(1:n, r$sesonal)

Проверка разделимости непериодических компонент + автогруппировка SSA

n <- 96*2-1
L <- 96
eps <- 1/97

C <- 1
omega_cs <- 1/12
omega_sn <- 1/24
a <- 1/100
f_sum <- function(x){
  f_const(x, C = C) +
    f_cos(x, omega = omega_cs) +
    f_exp(x, a = a) +
    f_sin(x, omega = omega_sn)
}


f_C <- f_const |> generate_ts(n, C = C)
f_c <- f_cos |> generate_ts(n, omega = omega_cs)
f_s <- f_sin |> generate_ts(n, omega = omega_sn)
f_e <- f_exp |> generate_ts(n, a = a)

f_n <- f_sum(1:n)

library(xtable)

Предупреждение: пакет ‘xtable’ был собран под R версии 4.2.3

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA", "CiSSA"),
  e_err = c(20, 20),
  c_err = c(23, 35),
  ec_err = c(20, 20),
  sin_err = c (20, 20),
  cos_err = c(1, 1)
)


# Отрисовка ряда f_n
plot(f_n, type = "l", lwd = 3, col = 'red', ylim = c(-2, 10),
     xlab = "Время", ylab = "Значения ряда", main = "Разложение временного ряда")

# Добавление отдельных компонентов (f_C, f_c, f_e)
lines(f_C, col = "blue")  # Компонент f_C

lines(f_c, col = "blue")  # Компонент f_c
lines(f_e, col = "blue")  # Компонент f_e

lines(f_s, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

c <- circulant_SSA(f_n, L = L, extend_flag = TRUE)
# r <- c$t_series
r <- grouping_cissa(c,
                    groups = list(
                      # trend = c(0, 1/100),
                      trend = c(-1, 1/25),
                      sesonal_cos = c(1/12-eps, 1/12+eps),
                      sesonal_sin = c(1/24 - eps, 1/24+eps)
                    ))$t_series

data$cos_err[2] <- mse(f_s, r$sesonal_sin) |> formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(f_c, r$sesonal_cos) |> formatC(format = "e", digits = 1)
data$ec_err[2] <- mse(f_C+f_e, r$trend) |> formatC(format = "e", digits = 1)


# png("C:/Users/nik1m/Desktop/уник/6 сем/курсач/Текст работы/img/trend inseparability/CiSSA.png")  # сохранение в формате PNG

plot(1:n, f_n, type = "l", lwd=3, ylim= c(-2, 10), col="red",
     xlab = "Время", ylab = "Значения ряда", main = "CiSSA разложение временного ряда")
lines(1:n, r$trend, col = "blue")

lines(1:n, r$sesonal_sin, col = "blue")
lines(1:n, r$sesonal_cos, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

# dev.off()  # завершение сохранения










s <- ssa(f_n, L)
#e <- eossa_new(s, nested.groups = list(1:10), clust_type = "distance", k= 7)
e <- eossa(s, 1:10, k = 7)

g_sesonal <- grouping.auto(e, base = "eigen",
                   freq.bins = list(trend = c(2*eps),
                                    sesonal2 = c(1/24-eps, 1/24+eps),
                                    sesonal1 = c(1/12-eps, 1/12+eps)
                                    ),
                   threshold = 0.1)


r <- Rssa::reconstruct(e, groups=c(list(exp = 1,
                                C = 2
                                ),
                             g_sesonal)
                 )

plot(wcor(e, groups = 1:24), scales = list(at = c(10, 20, 30)))

Предупреждение в wcor.ossa(e, groups = 1:24) :
  Component matrices are not F-orthogonal (max F-cor is 0.93). W-cor matrix can be irrelevant

data$c_err[1] <- mse(f_C, r$C) |> formatC(format = "e", digits = 1)
data$e_err[1] <- mse(f_e, r$exp) |> formatC(format = "e", digits = 1)
data$cos_err[1] <- mse(f_c, r$sesonal1) |> formatC(format = "e", digits = 1)
data$sin_err[1] <- mse(f_s, r$sesonal2) |> formatC(format = "e", digits = 1)
data$ec_err[1] <- mse(f_C+f_e, r$C+r$exp) |> formatC(format = "e", digits = 1)


# png("C:/Users/nik1m/Desktop/уник/6 сем/курсач/Текст работы/img/trend inseparability/SSA.png")  # сохранение в формате PNG

plot(1:n, f_n, type = "l", lwd=3, ylim= c(-2, 10), col="red",
     xlab = "Время", ylab = "Значения ряда", main = "SSA разложение временного ряда")

lines(1:n, r$trend, type = "l", col="green")

lines(1:n, r$exp, type = "l", ylim= c(-2, 10), col="blue")
lines(1:n, r$C, col = "blue")

lines(1:n, r$sesonal1, col = "blue")
lines(1:n, r$sesonal2, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

# Шаг 3: Преобразование данных в формат LaTeX
table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Fri Feb 28 16:59:21 2025
\begin{table}[ht]
\centering
\begin{tabular}{llllll}
  \hline
Метод & e\_err & c\_err & ec\_err & sin\_err & cos\_err \\ 
  \hline
SSA & 2.2e-25 & 2.2e-25 & 4.2e-28 & 3.8e-29 & 1.6e-29 \\ 
  CiSSA & 20 & 35 & 1.4e-03 & 7.7e-04 & 1.9e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Пример cos*exp

n <- 96*2-1
L <- 96
eps <- 1/(L+1)

C <- 1
omega_cs <- 1/12
omega_sn <- 1/24
a <- 1/100
omega_exp <- 1/48
f_sum <- function(x){
    f_cos(x, omega = omega_cs) +
    f_exp_mod_harm_series(x, a = a, omega = omega_exp) +
    f_sin(x, omega = omega_sn) 
}


f_c <- f_cos |> generate_ts(n, omega = omega_cs)
f_s <- f_sin |> generate_ts(n, omega = omega_sn)
f_e <- f_exp_mod_harm_series |> generate_ts(n, a = a, omega = omega_exp)

f_n <- f_sum(1:n)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA", "CiSSA"),
  exp_err = c(20, 20),
  sin_err = c (20, 20),
  cos_err = c(1, 1)
)


# Отрисовка ряда f_n
plot(f_n, type = "l", lwd = 3, col = 'red', ylim = c(-10, 10),
     xlab = "Время", ylab = "Значения ряда", main = "Разложение временного ряда")

# Добавление отдельных компонентов (f_C, f_c, f_e)
lines(f_c, col = "blue")  # Компонент f_c

lines(f_e, col = "blue")  # Компонент f_e
lines(f_s, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

c <- circulant_SSA(f_n, L = L, extend_flag = TRUE)
# r <- c$t_series
r <- grouping_cissa(c,
                    groups = list(
                      trend = c(0, 1/26-eps),
                      sesonal_cos = c(1/12-eps, 1/12+eps),
                      sesonal_sin = c(1/24-eps, 1/24+eps)
                    ))$t_series

data$cos_err[2] <- mse(f_s, r$sesonal_sin) |> formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(f_c, r$sesonal_cos) |> formatC(format = "e", digits = 1)
data$exp_err[2] <- mse(f_e, r$trend) |> formatC(format = "e", digits = 1)


# png("C:/Users/nik1m/Desktop/уник/6 сем/курсач/Текст работы/img/trend inseparability/CiSSA.png")  # сохранение в формате PNG

plot(1:n, f_n, type = "l", lwd=3, ylim= c(-10, 10), col="red",
     xlab = "Время", ylab = "Значения ряда", main = "CiSSA разложение временного ряда")
lines(1:n, r$trend, col = "blue")

lines(1:n, r$sesonal_sin, col = "blue")
lines(1:n, r$sesonal_cos, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

# dev.off()  # завершение сохранения










s <- ssa(f_n, L)
e <- eossa_new(s, nested.groups = list(1:30), clust_type = "distance")

g_sesonal <- grouping.auto(e, base = "eigen",
                   freq.bins = list(trend = c(1/25-eps),
                                    sesonal2 = c(1/24-eps, 1/24+eps),
                                    sesonal1 = c(1/12-eps, 1/12+eps)
                                    ),
                   threshold = 0.1)


r <- reconstruct(e, groups= g_sesonal)

plot(wcor(e, groups = 1:24), scales = list(at = c(10, 20, 30)))

Предупреждение в wcor.ossa(e, groups = 1:24) :
  Component matrices are not F-orthogonal (max F-cor is -0.529). W-cor matrix can be irrelevant

data$exp_err[1] <- mse(f_e, r$trend)  |> formatC(format = "e", digits = 1)
data$cos_err[1] <- mse(f_c, r$sesonal1) |> formatC(format = "e", digits = 1)
data$sin_err[1] <- mse(f_s, r$sesonal2) |> formatC(format = "e", digits = 1)


# png("C:/Users/nik1m/Desktop/уник/6 сем/курсач/Текст работы/img/trend inseparability/SSA.png")  # сохранение в формате PNG

plot(1:n, f_n, type = "l", lwd=3, ylim= c(-10, 10), col="red",
     xlab = "Время", ylab = "Значения ряда", main = "SSA разложение временного ряда")

lines(1:n, r$trend, type = "l", col="blue")

lines(1:n, r$sesonal1, col = "blue")
lines(1:n, r$sesonal2, col = "blue")

# Легенда
legend("topleft", legend = c("Весь ряд", "Компоненты"), 
       col = c("red", "blue"), lty = 1, lwd = 3)

# Шаг 3: Преобразование данных в формат LaTeX
table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Fri Feb 28 16:59:23 2025
\begin{table}[ht]
\centering
\begin{tabular}{llll}
  \hline
Метод & exp\_err & sin\_err & cos\_err \\ 
  \hline
SSA & 4.7e-29 & 1.1e-29 & 8.4e-30 \\ 
  CiSSA & 3.2e-02 & 1.0e-03 & 5.8e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Данные IP

library(readxl)

Предупреждение: пакет ‘readxl’ был собран под R версии 4.2.3

data <- read_excel("Data/International_Financial_Statistics_.xlsx")

New names:
• `` -> `...2`
• `` -> `...3`
• `` -> `...4`
• `` -> `...5`
• `` -> `...6`
• `` -> `...7`
• `` -> `...8`
• `` -> `...9`
• `` -> `...10`
• `` -> `...11`
• `` -> `...12`
• `` -> `...13`
• `` -> `...14`
• `` -> `...15`
• `` -> `...16`
• `` -> `...17`
• `` -> `...18`
• `` -> `...19`
• `` -> `...20`
• `` -> `...21`
• `` -> `...22`
• `` -> `...23`
• `` -> `...24`
• `` -> `...25`
• `` -> `...26`
• `` -> `...27`
• `` -> `...28`
• `` -> `...29`
• `` -> `...30`
• `` -> `...31`
• `` -> `...32`
• `` -> `...33`
• `` -> `...34`
• `` -> `...35`
• `` -> `...36`
• `` -> `...37`
• `` -> `...38`
• `` -> `...39`
• `` -> `...40`
• `` -> `...41`
• `` -> `...42`
• `` -> `...43`
• `` -> `...44`
• `` -> `...45`
• `` -> `...46`
• `` -> `...47`
• `` -> `...48`
• `` -> `...49`
• `` -> `...50`
• `` -> `...51`
• `` -> `...52`
• `` -> `...53`
• `` -> `...54`
• `` -> `...55`
• `` -> `...56`
• `` -> `...57`
• `` -> `...58`
• `` -> `...59`
• `` -> `...60`
• `` -> `...61`
• `` -> `...62`
• `` -> `...63`
• `` -> `...64`
• `` -> `...65`
• `` -> `...66`
• `` -> `...67`
• `` -> `...68`
• `` -> `...69`
• `` -> `...70`
• `` -> `...71`
• `` -> `...72`
• `` -> `...73`
• `` -> `...74`
• `` -> `...75`
• `` -> `...76`
• `` -> `...77`
• `` -> `...78`
• `` -> `...79`
• `` -> `...80`
• `` -> `...81`
• `` -> `...82`
• `` -> `...83`
• `` -> `...84`
• `` -> `...85`
• `` -> `...86`
• `` -> `...87`
• `` -> `...88`
• `` -> `...89`
• `` -> `...90`
• `` -> `...91`
• `` -> `...92`
• `` -> `...93`
• `` -> `...94`
• `` -> `...95`
• `` -> `...96`
• `` -> `...97`
• `` -> `...98`
• `` -> `...99`
• `` -> `...100`
• `` -> `...101`
• `` -> `...102`
• `` -> `...103`
• `` -> `...104`
• `` -> `...105`
• `` -> `...106`
• `` -> `...107`
• `` -> `...108`
• `` -> `...109`
• `` -> `...110`
• `` -> `...111`
• `` -> `...112`
• `` -> `...113`
• `` -> `...114`
• `` -> `...115`
• `` -> `...116`
• `` -> `...117`
• `` -> `...118`
• `` -> `...119`
• `` -> `...120`
• `` -> `...121`
• `` -> `...122`
• `` -> `...123`
• `` -> `...124`
• `` -> `...125`
• `` -> `...126`
• `` -> `...127`
• `` -> `...128`
• `` -> `...129`
• `` -> `...130`
• `` -> `...131`
• `` -> `...132`
• `` -> `...133`
• `` -> `...134`
• `` -> `...135`
• `` -> `...136`
• `` -> `...137`
• `` -> `...138`
• `` -> `...139`
• `` -> `...140`
• `` -> `...141`
• `` -> `...142`
• `` -> `...143`
• `` -> `...144`
• `` -> `...145`
• `` -> `...146`
• `` -> `...147`
• `` -> `...148`
• `` -> `...149`
• `` -> `...150`
• `` -> `...151`
• `` -> `...152`
• `` -> `...153`
• `` -> `...154`
• `` -> `...155`
• `` -> `...156`
• `` -> `...157`
• `` -> `...158`
• `` -> `...159`
• `` -> `...160`
• `` -> `...161`
• `` -> `...162`
• `` -> `...163`
• `` -> `...164`
• `` -> `...165`
• `` -> `...166`
• `` -> `...167`
• `` -> `...168`
• `` -> `...169`
• `` -> `...170`
• `` -> `...171`
• `` -> `...172`
• `` -> `...173`
• `` -> `...174`
• `` -> `...175`
• `` -> `...176`
• `` -> `...177`
• `` -> `...178`
• `` -> `...179`
• `` -> `...180`
• `` -> `...181`
• `` -> `...182`
• `` -> `...183`
• `` -> `...184`
• `` -> `...185`
• `` -> `...186`
• `` -> `...187`
• `` -> `...188`
• `` -> `...189`
• `` -> `...190`
• `` -> `...191`
• `` -> `...192`
• `` -> `...193`
• `` -> `...194`
• `` -> `...195`
• `` -> `...196`
• `` -> `...197`
• `` -> `...198`
• `` -> `...199`
• `` -> `...200`
• `` -> `...201`
• `` -> `...202`
• `` -> `...203`
• `` -> `...204`
• `` -> `...205`
• `` -> `...206`
• `` -> `...207`
• `` -> `...208`
• `` -> `...209`
• `` -> `...210`
• `` -> `...211`
• `` -> `...212`
• `` -> `...213`
• `` -> `...214`
• `` -> `...215`
• `` -> `...216`
• `` -> `...217`
• `` -> `...218`
• `` -> `...219`
• `` -> `...220`
• `` -> `...221`
• `` -> `...222`
• `` -> `...223`
• `` -> `...224`
• `` -> `...225`
• `` -> `...226`
• `` -> `...227`
• `` -> `...228`
• `` -> `...229`
• `` -> `...230`
• `` -> `...231`
• `` -> `...232`
• `` -> `...233`
• `` -> `...234`
• `` -> `...235`
• `` -> `...236`
• `` -> `...237`
• `` -> `...238`
• `` -> `...239`
• `` -> `...240`
• `` -> `...241`
• `` -> `...242`
• `` -> `...243`
• `` -> `...244`
• `` -> `...245`
• `` -> `...246`
• `` -> `...247`
• `` -> `...248`
• `` -> `...249`
• `` -> `...250`
• `` -> `...251`
• `` -> `...252`
• `` -> `...253`
• `` -> `...254`
• `` -> `...255`
• `` -> `...256`
• `` -> `...257`
• `` -> `...258`
• `` -> `...259`
• `` -> `...260`
• `` -> `...261`
• `` -> `...262`
• `` -> `...263`
• `` -> `...264`
• `` -> `...265`
• `` -> `...266`
• `` -> `...267`
• `` -> `...268`
• `` -> `...269`
• `` -> `...270`
• `` -> `...271`
• `` -> `...272`
• `` -> `...273`
• `` -> `...274`
• `` -> `...275`
• `` -> `...276`
• `` -> `...277`
• `` -> `...278`
• `` -> `...279`
• `` -> `...280`
• `` -> `...281`
• `` -> `...282`
• `` -> `...283`
• `` -> `...284`
• `` -> `...285`
• `` -> `...286`
• `` -> `...287`
• `` -> `...288`
• `` -> `...289`
• `` -> `...290`
• `` -> `...291`
• `` -> `...292`
• `` -> `...293`
• `` -> `...294`
• `` -> `...295`
• `` -> `...296`
• `` -> `...297`
• `` -> `...298`
• `` -> `...299`
• `` -> `...300`
• `` -> `...301`
• `` -> `...302`
• `` -> `...303`
• `` -> `...304`
• `` -> `...305`
• `` -> `...306`
• `` -> `...307`
• `` -> `...308`
• `` -> `...309`
• `` -> `...310`
• `` -> `...311`
• `` -> `...312`
• `` -> `...313`
• `` -> `...314`
• `` -> `...315`
• `` -> `...316`
• `` -> `...317`
• `` -> `...318`
• `` -> `...319`
• `` -> `...320`
• `` -> `...321`
• `` -> `...322`
• `` -> `...323`
• `` -> `...324`
• `` -> `...325`
• `` -> `...326`
• `` -> `...327`
• `` -> `...328`
• `` -> `...329`
• `` -> `...330`
• `` -> `...331`
• `` -> `...332`
• `` -> `...333`
• `` -> `...334`
• `` -> `...335`
• `` -> `...336`
• `` -> `...337`
• `` -> `...338`
• `` -> `...339`
• `` -> `...340`
• `` -> `...341`
• `` -> `...342`
• `` -> `...343`
• `` -> `...344`
• `` -> `...345`
• `` -> `...346`
• `` -> `...347`
• `` -> `...348`
• `` -> `...349`
• `` -> `...350`
• `` -> `...351`
• `` -> `...352`
• `` -> `...353`
• `` -> `...354`
• `` -> `...355`
• `` -> `...356`
• `` -> `...357`
• `` -> `...358`
• `` -> `...359`
• `` -> `...360`
• `` -> `...361`
• `` -> `...362`
• `` -> `...363`
• `` -> `...364`
• `` -> `...365`
• `` -> `...366`
• `` -> `...367`
• `` -> `...368`
• `` -> `...369`
• `` -> `...370`
• `` -> `...371`
• `` -> `...372`
• `` -> `...373`
• `` -> `...374`
• `` -> `...375`
• `` -> `...376`
• `` -> `...377`
• `` -> `...378`
• `` -> `...379`
• `` -> `...380`
• `` -> `...381`
• `` -> `...382`
• `` -> `...383`
• `` -> `...384`
• `` -> `...385`
• `` -> `...386`
• `` -> `...387`
• `` -> `...388`
• `` -> `...389`
• `` -> `...390`
• `` -> `...391`
• `` -> `...392`
• `` -> `...393`
• `` -> `...394`
• `` -> `...395`
• `` -> `...396`
• `` -> `...397`
• `` -> `...398`
• `` -> `...399`
• `` -> `...400`
• `` -> `...401`
• `` -> `...402`
• `` -> `...403`
• `` -> `...404`
• `` -> `...405`
• `` -> `...406`
• `` -> `...407`
• `` -> `...408`
• `` -> `...409`
• `` -> `...410`
• `` -> `...411`
• `` -> `...412`
• `` -> `...413`
• `` -> `...414`
• `` -> `...415`
• `` -> `...416`
• `` -> `...417`
• `` -> `...418`
• `` -> `...419`
• `` -> `...420`
• `` -> `...421`
• `` -> `...422`
• `` -> `...423`
• `` -> `...424`
• `` -> `...425`
• `` -> `...426`
• `` -> `...427`
• `` -> `...428`
• `` -> `...429`
• `` -> `...430`
• `` -> `...431`
• `` -> `...432`
• `` -> `...433`
• `` -> `...434`
• `` -> `...435`
• `` -> `...436`
• `` -> `...437`
• `` -> `...438`
• `` -> `...439`
• `` -> `...440`
• `` -> `...441`
• `` -> `...442`
• `` -> `...443`
• `` -> `...444`
• `` -> `...445`
• `` -> `...446`
• `` -> `...447`
• `` -> `...448`
• `` -> `...449`
• `` -> `...450`
• `` -> `...451`
• `` -> `...452`
• `` -> `...453`
• `` -> `...454`
• `` -> `...455`
• `` -> `...456`
• `` -> `...457`
• `` -> `...458`
• `` -> `...459`
• `` -> `...460`
• `` -> `...461`
• `` -> `...462`
• `` -> `...463`
• `` -> `...464`
• `` -> `...465`
• `` -> `...466`
• `` -> `...467`
• `` -> `...468`
• `` -> `...469`
• `` -> `...470`
• `` -> `...471`
• `` -> `...472`
• `` -> `...473`
• `` -> `...474`
• `` -> `...475`
• `` -> `...476`
• `` -> `...477`
• `` -> `...478`
• `` -> `...479`
• `` -> `...480`
• `` -> `...481`
• `` -> `...482`
• `` -> `...483`
• `` -> `...484`
• `` -> `...485`
• `` -> `...486`
• `` -> `...487`
• `` -> `...488`
• `` -> `...489`
• `` -> `...490`
• `` -> `...491`
• `` -> `...492`
• `` -> `...493`
• `` -> `...494`
• `` -> `...495`
• `` -> `...496`
• `` -> `...497`
• `` -> `...498`
• `` -> `...499`
• `` -> `...500`
• `` -> `...501`
• `` -> `...502`
• `` -> `...503`
• `` -> `...504`
• `` -> `...505`
• `` -> `...506`
• `` -> `...507`
• `` -> `...508`
• `` -> `...509`
• `` -> `...510`
• `` -> `...511`
• `` -> `...512`
• `` -> `...513`
• `` -> `...514`
• `` -> `...515`
• `` -> `...516`
• `` -> `...517`
• `` -> `...518`
• `` -> `...519`
• `` -> `...520`
• `` -> `...521`
• `` -> `...522`
• `` -> `...523`
• `` -> `...524`
• `` -> `...525`
• `` -> `...526`
• `` -> `...527`
• `` -> `...528`
• `` -> `...529`
• `` -> `...530`
• `` -> `...531`
• `` -> `...532`
• `` -> `...533`
• `` -> `...534`
• `` -> `...535`
• `` -> `...536`
• `` -> `...537`
• `` -> `...538`
• `` -> `...539`
• `` -> `...540`
• `` -> `...541`
• `` -> `...542`
• `` -> `...543`
• `` -> `...544`
• `` -> `...545`
• `` -> `...546`
• `` -> `...547`
• `` -> `...548`
• `` -> `...549`
• `` -> `...550`
• `` -> `...551`
• `` -> `...552`
• `` -> `...553`
• `` -> `...554`
• `` -> `...555`
• `` -> `...556`
• `` -> `...557`
• `` -> `...558`
• `` -> `...559`
• `` -> `...560`
• `` -> `...561`
• `` -> `...562`
• `` -> `...563`
• `` -> `...564`
• `` -> `...565`
• `` -> `...566`
• `` -> `...567`
• `` -> `...568`
• `` -> `...569`
• `` -> `...570`
• `` -> `...571`
• `` -> `...572`
• `` -> `...573`
• `` -> `...574`
• `` -> `...575`
• `` -> `...576`
• `` -> `...577`
• `` -> `...578`
• `` -> `...579`

data |> head()

Отрисовка данных IP

dates <- seq(as.Date("1970-01-01"), as.Date("2018-1-30"), by = "month")
IP_values <- data[2, -c(1, 2)] |> as.double() 
plot(dates, IP_values, type="l")

Cissa

Отрисовка трендовой составляющей чёрным цветом, основной временной ряд — красным

data_slice <- 1:537
dates_slice <- dates[data_slice]
IP_values_slice <- IP_values[data_slice]
eps <- 1/193

c <- circulant_SSA(IP_values_slice, L = 192, extend_flag = TRUE)
r <- c$t_series
r <- grouping_cissa(c,
                    groups = list(
                      trend = c(0, 1/192),
                      cycle = c(1/97, 5/95),
                      sesonal = c(1/13, 1/2+0.0001)
                    )
                    )$t_series
r_sesonal <-  grouping_cissa(c,
                             groups = list(
                              s1 = c(16/192 - eps, 16/192 + eps),
                              s2 = c(32/192 - eps, 32/192 + eps),
                              s3 = c(48/192 - eps, 48/192 + eps),
                              s4 = c(64/192 - eps, 64/192 + eps),
                              s5 = c(80/192 - eps, 80/192 + eps),
                              s6 = c(96/192 - eps, 96/192 + eps)
                             )
                             )$t_series
# cissa_trend <- r[,1] + r[,2]
# cissa_cycle <- r[, 3:11] |> rowSums()
# cissa_sesonal <- r[, c(17, 33, 49, 65, 81, 97)] |> rowSums()
# cissa_residuals <- IP_values_slice - (cissa_trend + cissa_cycle + cissa_sesonal)

cissa_trend <- r$trend
cissa_cycle <- r$cycle
cissa_sesonal <- Reduce("+", r_sesonal |> within(rm(residuals)))
cissa_residuals <- IP_values_slice - (cissa_trend + cissa_cycle + cissa_sesonal)


plot(dates_slice, IP_values_slice,
     type="l", col = "black")
lines(dates_slice, cissa_trend,
      type="l", col = "red")

plot(dates_slice, cissa_cycle,
     type="l", col = "red")

plot(dates_slice, cissa_sesonal,
     type="l", col = "red")

plot(dates_slice, cissa_residuals,
     type="l", col = "red")

plot(dates_slice, IP_values_slice,
     type="l", col = "black")
lines(dates_slice, cissa_trend+cissa_cycle+cissa_sesonal,
      type="l", col = "red")

SSA fossa

s <- ssa(IP_values_slice, L = 192)
e <- fossa(s)
# e <- eossa_new(s, nested.groups = list(1:30), clust_type = "distance")
eps <- 1/193

groups <- grouping.auto(e,
                   freq.bins = list(trend = c(1/192),
                                    cycle = c(1/97, 5/95),
                                    s1 = c(16/192 - eps, 16/192 + eps),
                                    s2 = c(32/192 - eps, 32/192 + eps),
                                    s3 = c(48/192 - eps, 48/192 + eps),
                                    s4 = c(64/192 - eps, 64/192 + eps),
                                    s5 = c(80/192 - eps, 80/192 + eps),
                                    s6 = c(96/192 - eps, 96/192 + eps)
                                    ),
                   threshold = 0)


plot(wcor(e, groups = 1:30), scales = list(at = c(10, 20, 30)),
     main = "W-correlation matrix SSA (fossa)")

r <- reconstruct(e, groups=groups)

ssa_trend_f <- r$trend
ssa_cycle_f <- r$cycle
ssa_sesonal_f <- r$s1 + r$s2 + r$s3 + r$s4 + r$s5 + r$s6
ssa_residuals_f <- IP_values_slice - (ssa_trend_f + ssa_cycle_f + ssa_sesonal_f)

plot(dates_slice, IP_values_slice,
     type="l", col = "black")
lines(dates_slice, ssa_trend_f,
      type="l", col = "magenta")

plot(dates_slice, ssa_cycle_f, 
     type="l", col = "magenta")

plot(dates_slice, ssa_sesonal_f, 
     type="l", col = "magenta")

plot(dates_slice, ssa_residuals_f,
     type="l", col = "magenta")

SSA eossa

library(Rssa)
source("eossa_new.r")
s <- ssa(IP_values_slice, L = 192)
e <- eossa_new(s, nested.groups = list(1:30), clust_type = "distance")




groups <- grouping.auto(e,
                   freq.bins = list(trend = c(1/192),
                                    cycle = c(1/97, 5/95),
                                    s1 = c(16/192 - eps, 16/192 + eps),
                                    s2 = c(32/192 - eps, 32/192 + eps),
                                    s3 = c(48/192 - eps, 48/192 + eps),
                                    s4 = c(64/192 - eps, 64/192 + eps),
                                    s5 = c(80/192 - eps, 80/192 + eps),
                                    s6 = c(96/192 - eps, 96/192 + eps)
                                    ),
                   threshold = 0)
plot(wcor(e, groups = 1:30), scales = list(at = c(10, 20, 30)),
     main = "W-correlation matrix SSA (eossa)")

Предупреждение в wcor.ossa(e, groups = 1:30) :
  Component matrices are not F-orthogonal (max F-cor is -0.0621). W-cor matrix can be irrelevant

r <- reconstruct(e, groups=groups)

ssa_trend <- r$trend
ssa_cycle <- r$cycle
ssa_sesonal <- r$s1 + r$s2 + r$s3 + r$s4 + r$s5 + r$s6
ssa_residuals <- IP_values_slice - (ssa_trend + ssa_cycle + ssa_sesonal)

plot(dates_slice, IP_values_slice,
     type="l", col = "black")
lines(dates_slice, ssa_trend,
      type="l", col = "blue")

plot(dates_slice, ssa_cycle, 
     type="l", col = "blue")

plot(dates_slice, ssa_sesonal, 
     type="l", col = "blue")

plot(dates_slice, ssa_residuals,
     type="l", col = "blue")

plot(dates_slice, IP_values_slice,
     main = "IP USA тренд",xlab = "Время", ylab = "Значение",
     type="l", col = "black")
lines(dates_slice, ssa_trend,
      type="l", col = "blue", lwd=2)

lines(dates_slice, ssa_trend_f,
      type="l", col = "magenta", lwd=2)
lines(dates_slice, cissa_trend,
      type="l", col = "red", lwd=2)

# Легенда
legend("topleft", legend = c("Весь ряд", "CiSSA тренд", "SSA тренд (eossa)", "SSA тренд (fossa)"), 
       col = c("black", "red", "blue", "magenta"), lty = 1, lwd = 3)

plot(dates_slice, ssa_cycle,
     main = "IP USA цикличность", xlab = "Время", ylab = "Значение",
     type="l", col = "blue", ylim=c(-10, 10), lwd=2)
lines(dates_slice, cissa_cycle,
      type="l", col = "red", lwd=2)

lines(dates_slice, ssa_cycle_f,
      type="l", col = "magenta", lwd=2)
# Легенда
legend("topleft", legend = c("CiSSA", "SSA (eossa)", "SSA (fossa)"), 
       col = c("red", "blue", "magenta"), lty = 1, lwd = 3)

plot(dates_slice, IP_values_slice,
     main = "IP USA тренд + цикличность",xlab = "Время", ylab = "Значение",
     type="l", col = "black")
lines(dates_slice, ssa_trend + ssa_cycle,
      type="l", col = "blue", lwd=2)

lines(dates_slice, ssa_trend_f + ssa_cycle_f,
      type="l", col = "magenta", lwd=2)
lines(dates_slice, cissa_trend + cissa_cycle,
      type="l", col = "red", lwd=2)

# Легенда
legend("topleft", legend = c("Весь ряд", "CiSSA", "SSA (eossa)", "SSA (fossa)"), 
       col = c("black", "red", "blue", "magenta"), lty = 1, lwd = 3)

NA
NA

# Настройка графиков для отображения двух графиков один под другим с общей осью X
layout(matrix(c(1, 2), nrow = 2, byrow = TRUE), heights = c(1, 1.2))

# Построение первого графика
par(mar = c(2, 4, 2, 2)) # Уменьшение нижнего отступа
plot(dates_slice, ssa_sesonal, type = "l", col = "blue", lwd = 1,
     main = "SSA (eossa) сезонность", xlab = "", ylab = "Значение")
# Добавление оси X внизу первого графика, но с пустыми метками
axis(1, labels = FALSE)

# Построение второго графика
par(mar = c(5, 4, 2, 2)) # Увеличение нижнего отступа
plot(dates_slice, ssa_sesonal_f, type = "l", col = "magenta", lwd = 1,
     main = "SSA (fossa) сезонность", xlab = "Время", ylab = "Значение")

par(mar = c(3, 4, 2, 2)) # Увеличение нижнего отступа

plot(dates_slice, cissa_sesonal, type = "l", col = "red", lwd = 1,
     main = "CiSSA сезонность", xlab = "Время", ylab = "Значение")

# Восстановление макета по умолчанию
layout(1)

NA
NA

plot(dates_slice, ssa_residuals, 
     main = "IP USA остаток", xlab = "Время", ylab = "Значение",
     type="l", col = "blue", ylim=c(-2, 2))
lines(dates_slice, cissa_residuals,
      type="l", col = "red")

lines(dates_slice, ssa_residuals_f,
      type="l", col = "magenta")
legend("topleft", legend = c("CiSSA", "SSA (eossa)", "SSA (fossa)"), 
       col = c("red", "blue", "magenta"), lty = 1, lwd = 3)

ssa_residuals |> density() |> plot()

cissa_residuals |> density() |> plot()

Отделение сигнала от шума

set.seed(100)

n_mse_tests <- function(n){
  n <- 96*2-1
  L <- 96
  sigma <- 0.1
  
  
  C <- 1
  omega_cs <- 1/12
  omega_sn <- 1/24
  a <- 1/100
  f_sum <- function(x){
    f_const(x, C = C) +
      f_cos(x, omega = omega_cs) +
      f_exp(x, a = a) +
      f_sin(x, omega = omega_sn)
  }
  
  
  f_C <- f_const |> generate_ts(n, C = C)
  f_c <- f_cos |> generate_ts(n, omega = omega_cs)
  f_s <- f_sin |> generate_ts(n, omega = omega_sn)
  f_e <- f_exp |> generate_ts(n, a = a)
  
  mse_lst <- list()
  for (i in 1:n) {
    f_noise <- rnorm(n, sd = sigma)
    
    f_n <- f_sum(1:n) + f_noise
    
    
    
    c <- circulant_SSA(f_n, L = L, extend_flag = TRUE)
    # r <- c$t_series
    r <- grouping_cissa(c, groups= list(trend = c(0, 1/1000), 
                                        sesonal2 = c(1/25, 1/23),
                                        sesonal1 = c(1/13, 1/10)
                                        ))$t_series
    
    # mse_lst$cissa <- c(mse_lst$cissa, mse(f_sum(1:n), r[, 9] + r[, 5] + r[, 1])) 
    mse_lst$cissa <- c(mse_lst$cissa,
                       mse(f_sum(1:n),
                           r$trend + r$sesonal1 + r$sesonal2)) 
    
    
    
    
    s <- ssa(f_n, L)
    # e <- eossa(s, 1:10, k = 6)
    e <- fossa(s)
    
    g_sesonal <- grouping.auto(e, base = "eigen",
                       freq.bins = list(trend = 1/1000, 
                                        sesonal2 = c(1/25, 1/23),
                                        sesonal1 = c(1/13, 1/10)
                                        ),
                       threshold = 0.5)
    
    r <- reconstruct(e, groups=c(list(exp = 1, C = 2), g_sesonal))
    
    mse_lst$ssa <- 
      c(mse_lst$ssa, mse(f_sum(1:n), r$trend + r$sesonal2 + r$sesonal1))
 
  }
  return(mse_lst)
}

res_mse_test <- n_mse_tests(100)

# Оценка плотности
density_estimate_cissa <- density(res_mse_test$cissa)

# Построение графика плотности
plot(density_estimate_cissa, main = "Оценка плотности", 
     xlab = "Значение", ylab = "Плотность", 
     col = "blue", lwd = 2)

density_estimate_ssa <- density(res_mse_test$ssa)

# Построение графика плотности
plot(density_estimate_ssa, main = "Оценка плотности", 
     xlab = "Значение", ylab = "Плотность", 
     col = "blue", lwd = 2)

res_mse_test$cissa |> summary()

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02299 0.03062 0.03345 0.03378 0.03575 0.04800 

res_mse_test$cissa |> sd()

[1] 0.004333975

res_mse_test$ssa |> summary()

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.000576 0.001710 0.002125 0.002228 0.002585 0.006311 

res_mse_test$ssa |> sd()

[1] 0.0008283841

Как выполняется расширение ряда

IP_values_slice |> extend(192) |>
  plot(type="l", lwd = 3, ylab="Значения", xlab = "Индексы", main="Расширение временного ряда IP values")
lines((192:(length(IP_values_slice) + 192 -1)) + 1, IP_values_slice, type="l", col="red")

n <- 96*2 + 6
x <- 0:(n-1)
y <- sin(2*pi/12 * x)
L <- 96+6
plot(y |> extend(L), type="l", lwd = 4, ylab="Значения", xlab = "Индексы", main="Расширение временного ряда sin")
lines((x+L+1), y, type="l", col="red", lwd = 2)

CiSSA через Фурье

IP_values

cissa_like_fourier_transform <- function(ts, L){
  reconstruct_fft <- function(x, y, frequencies) {
    fft_y <- fft(y)
  
    amplitudes <- Mod(fft_y)
    phases <- Arg(fft_y)
    
    reconstructed <- matrix(0, length(amplitudes), length(x))
    n <- length(amplitudes)
    for (i in 1:(length(amplitudes))) {
      reconstructed[i, ] <-
        amplitudes[i] * 
        cos(2 * pi * frequencies[i] * (x) + phases[i]) /
        n
    }
    return(list(
      ts = reconstructed,
      frequencies = frequencies
      ))
  }
  
  n <- length(ts)
  x <- 0:(n-1)
  L <- L
  frequencies <- (0:(L-1)) / L
  y <- ts
  y_main <- y
  x_main <- x
  
  X <- hankel(y, L)
  K <- dim(X)[2]
  res <- list()
  for (i in 1:K){
    y <- X[, i]
    x <- x_main[1:(1 + L - 1)]
    res[[i]] <- reconstruct_fft(x, y, frequencies)$ts
  }
  
  # for (i in 1:L){
  #   plot(x, res[[i]][9, ])
  # }
  
  res_mult <- res
  # res
  
  averaging <- function(res_comp_wise_mult){
    K <- dim(X)[2]
    counters <- rep(0, n)
    res <- matrix(0, ncol = n, nrow = L)
    for (i in 1:length(res_comp_wise_mult)){
      res[, i:(i+L-1)] <- res[, i:(i+L-1)] + res_comp_wise_mult[[i]]
      counters[i:(i+L-1)] <- counters[i:(i+L-1)] + 1
    }
    for (i in 1:n){
      res[, i] <- res[, i] / counters[i]
    }
    res
  }
  
  avr <- averaging(res_mult)
  
  group_by_elementary_freq_foureir <- function(res_averaged){
    nf2 <- 0
    if (L %% 2) {
      nf2 <- (L + 1) / 2 - 1
    } else {
      nf2 <- L / 2 - 1
    }
    nft <- nf2 + abs((L %% 2) - 2)
    
    Z <- matrix(0, ncol = nft, nrow = n)
    
    # print(Z |> dim())
    # print(res_averaged |> dim())
    
    Z[, 1] <- res_averaged[1, ]
    for (k in 1:nf2) {
      Z[, k + 1] <- res_averaged[k + 1, ] + res_averaged[L + 2 - (k + 1), ]
    }
    if (L %% 2 != 0) {
      Z[, nft] <- res_averaged[nft, ]
    }
    
    
    return(list(
      t_series = Z,
      freq = (0:dim(Z)[2])/L
    ))
  }


  rs <- group_by_elementary_freq_foureir(avr)
  return(rs)
}

data_slice <- 1:537
dates_slice <- dates[data_slice]
IP_values_slice <- IP_values[data_slice]
eps <- 1/193


c <- circulant_SSA(IP_values_slice, L = 192, extend_flag = FALSE)
r <- c$t_series

c_ft <- cissa_like_fourier_transform(IP_values_slice, L = 192)
r_ft <- c_ft$t_series

for (i in 1:dim(r)[2]){
  plot(data_slice, r_ft[, i], col= "red", type = "l", lwd = 2)
  lines(data_slice, r[, i])
}

Сравнение разложений всего временного ряда через метод Фурье, SSA, auto_SSA, CiSSA и CiSSA с расширением ряда

reconstruct_fft <- function(x_init, y_init, extend_flag = FALSE) {
    x <- x_init
    y <- y_init
    N <- length(y_init)
    H <- 0
    if (extend_flag == TRUE){
      H <- length(y) %/% 2
      y <- extend(y, H)
      x <- 0:(length(y) - 1)
    }
    
    frequencies <- (0:(length(x)-1)) / length(x)
    fft_y <- fft(y)
    
    amplitudes <- Mod(fft_y)
    phases <- Arg(fft_y)
    
    reconstructed <- matrix(0, length(amplitudes), length(x))
    n <- length(amplitudes)
    L <- n
    for (i in 1:(length(amplitudes))) {
      reconstructed[i, ] <-
        amplitudes[i] * 
        cos(2 * pi * frequencies[i] * (x) + phases[i]) /
        n
    }
    
    
    
  
    group_by_elementary_freq_foureir <- function(res_averaged){
      nf2 <- 0
      if (L %% 2) {
        nf2 <- (L + 1) / 2 - 1
      } else {
        nf2 <- L / 2 - 1
      }
      nft <- nf2 + abs((L %% 2) - 2)
      
      Z <- matrix(0, ncol = nft, nrow = n)
      
      # print(Z |> dim())
      # print(res_averaged |> dim())
      
      Z[, 1] <- res_averaged[1, ]
      for (k in 1:nf2) {
        Z[, k + 1] <- res_averaged[k + 1, ] + res_averaged[L + 2 - (k + 1), ]
      }
      if (L %% 2 != 0) {
        Z[, nft] <- res_averaged[nft, ]
      }
      
      
      return(list(
        t_series = Z[(H+1):(N+H), ],
        freq = (0:(dim(Z)[2]-1))/L
      ))
    }
  
  
    rs <- group_by_elementary_freq_foureir(reconstructed)
    
    
    return(list(
      t_series = rs$t_series,
      freq = rs$freq
      ))
  }

Идеальный случай для CiSSA и Fourier

n <- 96*2
L <- 96
x <- 0:(n-1)
y1 <- sin(2*pi/12 * x)
y2 <- cos(2*pi/3 * x)/2
y <- y1 + y2
X <- hankel(y, L = L)
eps <- 1/(n+1)

s_ssa <- ssa(y[1:(n-1)], L)
r_ssa <- reconstruct(s_ssa, groups = list(
  F1 = c(1, 2),
  F2 = c(3, 4)
))
# r_ssa <- reconstruct(s_ssa, groups=list(
#   sesonal_sin = c(1, 2),
#   sesonal_cos = c(3, 4)
# ))
# plot(x, r_ssa$F1, type="l")
# plot(x, r_ssa$F2, type= "l")
# plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
e_ssa <- eossa_new(s_ssa, nested.groups = list(1:4), clust_type = "distance")
g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                       freq.bins = list(
                         sesonal_sin = c(1/12-eps, 1/12+eps),
                         sesonal_cos = c(1/3-eps, 1/3+eps)
                                        ),
                       threshold = 0.5)
r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
# plot(x, r_ssa_e$sesonal_sin)
# plot(x, r_ssa_e$sesonal_cos)



r_fft <- reconstruct_fft(x, y)
r_fft_grouped <- grouping_cissa(r_fft,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

r_fft_extended <- reconstruct_fft(x, y, TRUE)
r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# r_fft_grouped$sesonal_sin |> length()

# plot(x, r_fft_grouped$sesonal_sin)
# plot(x, r_fft_grouped$sesonal_cos)
# plot(x, r_fft_grouped$residuals)


r_cissa <- circulant_SSA(y, L)
r_cissa_grouped <- grouping_cissa(r_cissa,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series


r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# plot(x, r_cissa_grouped$sesonal_sin)
# plot(x, r_cissa_grouped$sesonal_cos)
# plot(x, r_cissa_grouped$residuals)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = c (1, 20, 20, 1, 1, 1),
  cos_err = c(1, 1, 1, 1, 1, 1),
  all_err = c(1, 1, 1, 1, 1, 1)
)

data$sin_err[1] <- mse(y1[1:(n-1)], r_ssa$F1) |>
  formatC(format = "e", digits = 1) 
data$cos_err[1] <- mse(y2[1:(n-1)], r_ssa$F2) |>
  formatC(format = "e", digits = 1)
data$all_err[1] <- mse(y1[1:(n-1)] + y2[1:(n-1)], r_ssa$F1 + r_ssa$F2) |>
  formatC(format = "e", digits = 1)

data$cos_err[2] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$all_err[2] <- mse(y1[1:(n-1)] + y2[1:(n-1)],
                       r_ssa_e$sesonal_sin + r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mse(y1, r_fft_grouped$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mse(y2, r_fft_grouped$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$all_err[3] <- mse(y1 + y2,
                       r_fft_grouped$sesonal_sin + r_fft_grouped$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mse(y1, r_cissa_grouped$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mse(y2, r_cissa_grouped$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$all_err[4] <- mse(y1 + y2,
                       r_cissa_grouped$sesonal_sin + r_cissa_grouped$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[5] <- mse(y1, r_cissa_grouped_ext$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mse(y2, r_cissa_grouped_ext$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$all_err[5] <- mse(y1 + y2,
                       r_cissa_grouped_ext$sesonal_sin + r_cissa_grouped_ext$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[6] <- mse(y1, r_fft_grouped_extended$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mse(y2, r_fft_grouped_extended$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$all_err[6] <- mse(y1 + y2,
                       r_fft_grouped_extended$sesonal_sin + r_fft_grouped_extended$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)


table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Fri Feb 28 17:00:57 2025
\begin{table}[ht]
\centering
\begin{tabular}{llll}
  \hline
Метод & sin\_err & cos\_err & all\_err \\ 
  \hline
SSA & 6.8e-30 & 1.5e-29 & 1.8e-29 \\ 
  SSA EOSSA & 1.5e-29 & 7.5e-30 & 2.0e-29 \\ 
  Fourier & 1.7e-28 & 3.5e-28 & 5.1e-28 \\ 
  CiSSA & 1.9e-29 & 5.3e-30 & 2.1e-29 \\ 
  CiSSA extended & 2.0e-04 & 8.6e-04 & 1.1e-03 \\ 
  Fourier extended & 6.2e-04 & 2.6e-03 & 3.2e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Идеальный случай с шумом

data <- list(
  Метод = c("SSA","SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = list(list(), list(), list(), list(), list(), list()),
  cos_err = list(list(), list(), list(), list(), list(), list()),
  all_err = list(list(), list(), list(), list(), list(), list())
)

for (i in 1:100){
  set.seed(i)
  
  n <- 96*2
  L <- 96
  x <- 0:(n-1)
  y1 <- sin(2*pi/12 * x)
  y2 <- cos(2*pi/3 * x)/2
  y <- y1 + y2 + rnorm(n, 0, 0.1)
  X <- hankel(y, L = L)
  eps <- 1/(n+1)
  
  s_ssa <- ssa(y[1:(n-1)], L)
  r_ssa <- reconstruct(s_ssa, groups = list(
    F1 = c(1, 2),
    F2 = c(3, 4)
  ))
  # r_ssa <- reconstruct(s_ssa, groups=list(
  #   sesonal_sin = c(1, 2),
  #   sesonal_cos = c(3, 4)
  # ))
  # plot(x[1:(n-1)], r_ssa$F1, type = "l")
  # plot(x[1:(n-1)], r_ssa$F2, type = "l")
  # plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
  e_ssa <- eossa_new(s_ssa, nested.groups = list(1:4), clust_type = "distance")
  g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                         freq.bins = list(
                           sesonal_sin = c(1/12-eps, 1/12+eps),
                           sesonal_cos = c(1/3-eps, 1/3+eps)
                                          ),
                         threshold = 0.5)
  r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
  # plot(x, r_ssa_e$sesonal_sin)
  # plot(x, r_ssa_e$sesonal_cos)
  
  
  
  r_fft <- reconstruct_fft(x, y)
  r_fft_grouped <- grouping_cissa(r_fft,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  r_fft_extended <- reconstruct_fft(x, y, TRUE)
  r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  # r_fft_grouped$sesonal_sin |> length()
  
  # plot(x, r_fft_grouped$sesonal_sin)
  # plot(x, r_fft_grouped$sesonal_cos)
  # plot(x, r_fft_grouped$residuals)
  
  
  r_cissa <- circulant_SSA(y, L)
  r_cissa_grouped <- grouping_cissa(r_cissa,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  
  r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
  r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  data$cos_err[[1]][[i]] <- 
    min(mse(y1[1:(n-1)], r_ssa$F1), mse(y1[1:(n-1)], r_ssa$F2))
  data$sin_err[[1]][[i]] <-
    min(mse(y2[1:(n-1)], r_ssa$F1), mse(y2[1:(n-1)], r_ssa$F2))
  data$all_err[[1]][[i]] <-
    min(mse(y1[1:(n-1)]+y2[1:(n-1)], r_ssa$F1+ r_ssa$F2))
  
  data$cos_err[[2]][[i]] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin)
  data$sin_err[[2]][[i]] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos)
  data$all_err[[2]][[i]] <- mse(y1[1:(n-1)]+ y2[1:(n-1)],
                                r_ssa_e$sesonal_sin + r_ssa_e$sesonal_cos)
  
  
  data$cos_err[[3]][[i]] <- mse(y1, r_fft_grouped$sesonal_sin)
  data$sin_err[[3]][[i]] <- mse(y2, r_fft_grouped$sesonal_cos)
  data$all_err[[3]][[i]] <- mse(y1 + y2,
                                r_fft_grouped$sesonal_cos +  r_fft_grouped$sesonal_sin)
  
  
  data$cos_err[[4]][[i]] <- mse(y1, r_cissa_grouped$sesonal_sin)
  data$sin_err[[4]][[i]] <- mse(y2, r_cissa_grouped$sesonal_cos)
  data$all_err[[4]][[i]] <- mse(y1 + y2,
                          r_cissa_grouped$sesonal_sin +  r_cissa_grouped$sesonal_cos
                                )
  
  data$cos_err[[5]][[i]] <- mse(y1, r_cissa_grouped_ext$sesonal_sin)
  data$sin_err[[5]][[i]] <- mse(y2, r_cissa_grouped_ext$sesonal_cos)
  data$all_err[[5]][[i]] <- mse(y1 + y2,
                          r_cissa_grouped_ext$sesonal_sin + r_cissa_grouped_ext$sesonal_cos
                                )
  
  data$cos_err[[6]][[i]] <- mse(y1, r_fft_grouped_extended$sesonal_sin)
  data$sin_err[[6]][[i]] <- mse(y2, r_fft_grouped_extended$sesonal_cos)
  data$all_err[[6]][[i]] <- mse(y1 + y2,
                          r_fft_grouped_extended$sesonal_sin + r_fft_grouped_extended$sesonal_cos
                                )
}
library(xtable)

# Шаг 2: Создание примера данных
data_prev <- data

data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = c(0, 0, 0, 0, 0, 0),
  cos_err = c(0, 0, 0, 0, 0, 0),
  all_err = c(0, 0, 0, 0, 0, 0)
)

data$cos_err[1] <- mean(data_prev$cos_err[[1]] |> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mean(data_prev$sin_err[[1]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[1] <- mean(data_prev$all_err[[1]]|> unlist()) |>
  formatC(format = "e", digits = 1)

data$cos_err[2] <- mean(data_prev$cos_err[[2]] |> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mean(data_prev$sin_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[2] <- mean(data_prev$all_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mean(data_prev$cos_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mean(data_prev$sin_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[3] <- mean(data_prev$all_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mean(data_prev$cos_err[[4]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mean(data_prev$sin_err[[4]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$all_err[4] <- mean(data_prev$all_err[[4]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[5] <- mean(data_prev$cos_err[[5]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mean(data_prev$sin_err[[5]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$all_err[5] <- mean(data_prev$all_err[[5]]|> unlist()) |>
  formatC(format = "e", digits = 1)

data$cos_err[6] <- mean(data_prev$cos_err[[6]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mean(data_prev$sin_err[[6]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$all_err[6] <- mean(data_prev$all_err[[6]]|> unlist()) |>
  formatC(format = "e", digits = 1)

for (i in 1:6){
  for (j in (i):6){
    x <- data_prev$cos_err[[i]] |> unlist()
    y <- data_prev$cos_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    print(paste("cos, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))
  }
}

[1] "cos,  SSA   SSA , p-val =  NaN"
[1] "cos,  SSA   SSA EOSSA , p-val =  4.4705711361447e-05"
[1] "cos,  SSA   Fourier , p-val =  5.30848050504907e-20"
[1] "cos,  SSA   CiSSA , p-val =  2.98626448986158e-20"
[1] "cos,  SSA   CiSSA extended , p-val =  5.35308949018552e-37"
[1] "cos,  SSA   Fourier extended , p-val =  2.61456030462997e-59"
[1] "cos,  SSA EOSSA   SSA EOSSA , p-val =  NaN"
[1] "cos,  SSA EOSSA   Fourier , p-val =  5.30815210685059e-20"
[1] "cos,  SSA EOSSA   CiSSA , p-val =  2.98721789488757e-20"
[1] "cos,  SSA EOSSA   CiSSA extended , p-val =  5.24897277920756e-37"
[1] "cos,  SSA EOSSA   Fourier extended , p-val =  2.58891311835391e-59"
[1] "cos,  Fourier   Fourier , p-val =  NaN"
[1] "cos,  Fourier   CiSSA , p-val =  1.70228094341221e-10"
[1] "cos,  Fourier   CiSSA extended , p-val =  4.47998903141858e-42"
[1] "cos,  Fourier   Fourier extended , p-val =  2.90661284766219e-62"
[1] "cos,  CiSSA   CiSSA , p-val =  NaN"
[1] "cos,  CiSSA   CiSSA extended , p-val =  3.33436422220891e-40"
[1] "cos,  CiSSA   Fourier extended , p-val =  3.27079203462158e-61"
[1] "cos,  CiSSA extended   CiSSA extended , p-val =  NaN"
[1] "cos,  CiSSA extended   Fourier extended , p-val =  1.82605722656699e-64"
[1] "cos,  Fourier extended   Fourier extended , p-val =  NaN"

for (i in 1:6){
  for (j in (i):6){
    x <- data_prev$sin_err[[i]] |> unlist()
    y <- data_prev$sin_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    print(paste("sin, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))
  }
}

[1] "sin,  SSA   SSA , p-val =  NaN"
[1] "sin,  SSA   SSA EOSSA , p-val =  5.55731118903082e-18"
[1] "sin,  SSA   Fourier , p-val =  7.32586799185552e-20"
[1] "sin,  SSA   CiSSA , p-val =  5.06061567708895e-19"
[1] "sin,  SSA   CiSSA extended , p-val =  1.47978944226882e-21"
[1] "sin,  SSA   Fourier extended , p-val =  7.90332165341783e-43"
[1] "sin,  SSA EOSSA   SSA EOSSA , p-val =  NaN"
[1] "sin,  SSA EOSSA   Fourier , p-val =  7.00903680548176e-20"
[1] "sin,  SSA EOSSA   CiSSA , p-val =  4.86142494601895e-19"
[1] "sin,  SSA EOSSA   CiSSA extended , p-val =  1.75054970045824e-21"
[1] "sin,  SSA EOSSA   Fourier extended , p-val =  8.8573857160748e-43"
[1] "sin,  Fourier   Fourier , p-val =  NaN"
[1] "sin,  Fourier   CiSSA , p-val =  1.84578267967756e-11"
[1] "sin,  Fourier   CiSSA extended , p-val =  3.59520276363416e-35"
[1] "sin,  Fourier   Fourier extended , p-val =  3.2131198725477e-51"
[1] "sin,  CiSSA   CiSSA , p-val =  NaN"
[1] "sin,  CiSSA   CiSSA extended , p-val =  3.86571664691894e-32"
[1] "sin,  CiSSA   Fourier extended , p-val =  4.37350616684091e-50"
[1] "sin,  CiSSA extended   CiSSA extended , p-val =  NaN"
[1] "sin,  CiSSA extended   Fourier extended , p-val =  4.6461440598169e-38"
[1] "sin,  Fourier extended   Fourier extended , p-val =  NaN"

for (i in 1:6){
  for (j in (i):6){
    x <- data_prev$all_err[[i]] |> unlist()
    y <- data_prev$all_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    print(paste("all_err, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))
  }
}

[1] "all_err,  SSA   SSA , p-val =  NaN"
[1] "all_err,  SSA   SSA EOSSA , p-val =  0.217780725368394"
[1] "all_err,  SSA   Fourier , p-val =  2.297959821497e-29"
[1] "all_err,  SSA   CiSSA , p-val =  1.79607090075285e-31"
[1] "all_err,  SSA   CiSSA extended , p-val =  2.17727486455347e-41"
[1] "all_err,  SSA   Fourier extended , p-val =  1.59940501680576e-65"
[1] "all_err,  SSA EOSSA   SSA EOSSA , p-val =  NaN"
[1] "all_err,  SSA EOSSA   Fourier , p-val =  2.29795982149755e-29"
[1] "all_err,  SSA EOSSA   CiSSA , p-val =  1.796070900754e-31"
[1] "all_err,  SSA EOSSA   CiSSA extended , p-val =  2.17727486455359e-41"
[1] "all_err,  SSA EOSSA   Fourier extended , p-val =  1.5994050168058e-65"
[1] "all_err,  Fourier   Fourier , p-val =  NaN"
[1] "all_err,  Fourier   CiSSA , p-val =  5.10461358074915e-16"
[1] "all_err,  Fourier   CiSSA extended , p-val =  3.96748674539869e-50"
[1] "all_err,  Fourier   Fourier extended , p-val =  6.13502561191831e-70"
[1] "all_err,  CiSSA   CiSSA , p-val =  NaN"
[1] "all_err,  CiSSA   CiSSA extended , p-val =  8.73138603765015e-48"
[1] "all_err,  CiSSA   Fourier extended , p-val =  5.59699676025336e-69"
[1] "all_err,  CiSSA extended   CiSSA extended , p-val =  NaN"
[1] "all_err,  CiSSA extended   Fourier extended , p-val =  1.1567123543996e-69"
[1] "all_err,  Fourier extended   Fourier extended , p-val =  NaN"

table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Fri Feb 28 17:01:49 2025
\begin{table}[ht]
\centering
\begin{tabular}{llll}
  \hline
Метод & sin\_err & cos\_err & all\_err \\ 
  \hline
SSA & 2.9e-04 & 3.1e-04 & 5.9e-04 \\ 
  SSA EOSSA & 2.9e-04 & 3.1e-04 & 5.9e-04 \\ 
  Fourier & 1.0e-04 & 1.1e-04 & 2.2e-04 \\ 
  CiSSA & 1.6e-04 & 1.8e-04 & 3.4e-04 \\ 
  CiSSA extended & 6.6e-04 & 1.9e-03 & 2.5e-03 \\ 
  Fourier extended & 1.2e-03 & 3.9e-03 & 5.1e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

# mean(data_prev$sin_err[[2]] |> unlist()) |> print()
# mean(data_prev$sin_err[[1]] |> unlist()) |> print()

Добавим непериодичность

n <- 96*2
x <- 0:(n-1)
L <- 96
y1 <- sin(2*pi/12 * x)
y2 <- cos(2*pi/3 * x)/2
y3 <- exp(x/100) + 1
y <- y1 + y2 + y3
eps <- 1/(n+1)

s_ssa <- ssa(y[1:(n-1)], L)
r_ssa <- reconstruct(s_ssa, groups=list(
  e = c(1),
  sesonal_sin = c(2, 3),
  sesonal_cos = c(4, 5)
))
# plot(x[1:(n-1)], r_ssa$e, type = "l")
# plot(x[1:(n-1)], r_ssa$sesonal_sin, type = "l")
# plot(x[1:(n-1)], r_ssa$sesonal_cos, type = "l")
# plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
e_ssa <- eossa_new(s_ssa, nested.groups = list(1:7), clust_type = "distance")
g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                       freq.bins = list(
                         e = c(0, 1/12-eps-eps),
                         sesonal_sin = c(1/12-eps, 1/12+eps),
                         sesonal_cos = c(1/3-eps, 1/3+eps)
                                        ),
                       threshold = 0.5)
r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
# plot(x, r_ssa_e$sesonal_sin)
# plot(x, r_ssa_e$sesonal_cos)



r_fft <- reconstruct_fft(x, y)
r_fft_grouped <- grouping_cissa(r_fft,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

r_fft_extended <- reconstruct_fft(x, y, TRUE)
r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# r_fft_grouped$sesonal_sin |> length()

# plot(x, r_fft_grouped$sesonal_sin)
# plot(x, r_fft_grouped$sesonal_cos)
# plot(x, r_fft_grouped$residuals)


r_cissa <- circulant_SSA(y, L)
r_cissa_grouped <- grouping_cissa(r_cissa,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series


r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# plot(x, r_cissa_grouped$sesonal_sin)
# plot(x, r_cissa_grouped$sesonal_cos)
# plot(x, r_cissa_grouped$residuals)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  e_err = c(1, 1, 1, 1, 1, 1),
  sin_err = c (1, 20, 20, 1, 1, 1),
  cos_err = c(1, 1, 1, 1, 1, 1),
  all_err = c(1, 1, 1, 1, 1, 1)
)

data$cos_err[1] <- mse(y1[1:(n-1)], r_ssa$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mse(y2[1:(n-1)], r_ssa$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[1] <- mse(y3[1:(n-1)], r_ssa$e) |>
  formatC(format = "e", digits = 1)
data$all_err[1] <- mse(y1[1:(n-1)]+ y2[1:(n-1)] + y3[1:(n-1)],
                     r_ssa$sesonal_sin + r_ssa$sesonal_cos + r_ssa$e) |>
  formatC(format = "e", digits = 1)

data$cos_err[2] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[2] <- mse(y3[1:(n-1)], r_ssa_e$e) |>
  formatC(format = "e", digits = 1)
data$all_err[2] <- mse(y1[1:(n-1)]+ y2[1:(n-1)] + y3[1:(n-1)],
                     r_ssa_e$sesonal_sin + r_ssa_e$sesonal_cos + r_ssa_e$e
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mse(y1, r_fft_grouped$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mse(y2, r_fft_grouped$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[3] <- mse(y3, r_fft_grouped$e) |>
  formatC(format = "e", digits = 1)
data$all_err[3] <- mse(y1 + y2 + y3, 
                       r_fft_grouped$sesonal_sin + r_fft_grouped$sesonal_cos + r_fft_grouped$e) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mse(y1, r_cissa_grouped$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mse(y2, r_cissa_grouped$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[4] <- mse(y3, r_cissa_grouped$e) |>
  formatC(format = "e", digits = 1)
data$all_err[4] <- mse(y1 + y2 + y3, 
                       r_cissa_grouped$sesonal_sin + 
                         r_cissa_grouped$sesonal_cos + 
                         r_cissa_grouped$e
                       ) |>
  formatC(format = "e", digits = 1)

data$cos_err[5] <- mse(y1, r_cissa_grouped_ext$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mse(y2, r_cissa_grouped_ext$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[5] <- mse(y3, r_cissa_grouped_ext$e) |>
  formatC(format = "e", digits = 1)
data$all_err[5] <- mse(y1 + y2 + y3, 
                       r_cissa_grouped_ext$sesonal_sin + 
                         r_cissa_grouped_ext$sesonal_cos +
                          r_cissa_grouped_ext$e
                       ) |>
  formatC(format = "e", digits = 1)

data$cos_err[6] <- mse(y1, r_fft_grouped_extended$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mse(y2, r_fft_grouped_extended$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[6] <- mse(y3, r_fft_grouped_extended$e) |>
  formatC(format = "e", digits = 1)
data$all_err[6] <- mse(y1 + y2 + y3, 
                       r_fft_grouped_extended$sesonal_sin +
                         r_fft_grouped_extended$sesonal_cos +
                         r_fft_grouped_extended$e
                       ) |>
  formatC(format = "e", digits = 1)



table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Fri Feb 28 17:01:50 2025
\begin{table}[ht]
\centering
\begin{tabular}{lllll}
  \hline
Метод & e\_err & sin\_err & cos\_err & all\_err \\ 
  \hline
SSA & 5.0e-03 & 8.9e-07 & 5.2e-05 & 4.4e-03 \\ 
  SSA EOSSA & 1.7e-28 & 1.6e-29 & 8.7e-30 & 1.6e-28 \\ 
  Fourier & 1.1e-01 & 6.1e-04 & 6.8e-03 & 1.1e-01 \\ 
  CiSSA & 5.3e-02 & 1.6e-05 & 4.9e-04 & 4.4e-02 \\ 
  CiSSA extended & 5.0e-04 & 2.1e-04 & 1.1e-03 & 6.0e-04 \\ 
  Fourier extended & 1.4e-03 & 1.3e-03 & 8.4e-03 & 9.6e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Непериодичность с шумом

# data <- data.frame(
#   Метод = c("Fourier", "CiSSA", "CiSSA с расширением ряда"),
#   sin_err = c (20, 20, 20),
#   cos_err = c(1, 1, 20),
#   exp_err = c(1, 1, 20)
# )
data <- list(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = list(list(), list(), list(), list(), list(), list()),
  cos_err = list(list(), list(), list(), list(), list(), list()),
  exp_err = list(list(), list(), list(), list(), list(), list()),
  all_err = list(list(), list(), list(), list(), list(), list())
)

for (i in 1:100){
  set.seed(i)
  
  n <- 96*2
  x <- 0:(n-1)
  L <- 96
  y1 <- sin(2*pi/12 * x)
  y2 <- cos(2*pi/3 * x)/2
  y3 <- exp(x/100) + 1
  y <- y1 + y2 + y3 + rnorm(n, 0, 0.1)
  eps <- 1/(n+1)
  
  s_ssa <- ssa(y[1:(n-1)], L)
  r_ssa <- reconstruct(s_ssa, groups=list(
    e = 1,
    sesonal_sin = c(2, 3),
    sesonal_cos = c(4, 5)
  ))
  
  
  e_ssa <- eossa_new(s_ssa, nested.groups = list(1:7), clust_type = "distance")
  g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                         freq.bins = list(
                           e = c(0, 1/12-eps-eps),
                           sesonal_sin = c(1/12-eps, 1/12+eps),
                           sesonal_cos = c(1/3-eps, 1/3+eps)
                                          ),
                         threshold = 0.5)
  r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
  
  
  
  r_fft <- reconstruct_fft(x, y)
  r_fft_grouped <- grouping_cissa(r_fft,
                      groups = list(
                        e = c(0, 1/12-eps-eps),
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  r_fft_extended <- reconstruct_fft(x, y, TRUE)
  r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                      groups = list(
                        e = c(0, 1/12-eps-eps),
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  
  r_cissa <- circulant_SSA(y, L)
  r_cissa_grouped <- grouping_cissa(r_cissa,
                      groups = list(
                        e = c(0, 1/12-eps-eps),
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  
  r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
  r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                      groups = list(
                        e = c(0, 1/12-eps-eps),
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  
  
  data$cos_err[[1]][[i]] <- mse(y1[1:(n-1)], r_ssa$sesonal_sin)
  data$sin_err[[1]][[i]] <- mse(y2[1:(n-1)], r_ssa$sesonal_cos)
  data$exp_err[[1]][[i]] <- mse(y3[1:(n-1)], r_ssa$e)
  data$all_err[[1]][[i]] <- mse(y1[1:(n-1)] + y2[1:(n-1)] + y3[1:(n-1)],
                                r_ssa$sesonal_sin+
                                  r_ssa$sesonal_cos +
                                  r_ssa$e
                                )
  
  data$cos_err[[2]][[i]] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin)
  data$sin_err[[2]][[i]] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos)
  data$exp_err[[2]][[i]] <- mse(y3[1:(n-1)], r_ssa_e$e)
  data$all_err[[2]][[i]] <- mse(y1[1:(n-1)] + y2[1:(n-1)] + y3[1:(n-1)],
                                r_ssa_e$sesonal_sin+ 
                                  r_ssa_e$sesonal_cos + 
                                  r_ssa_e$e
                                )
  
  
  # print(data$sin_err[[2]][[i]])
  
  
  data$cos_err[[3]][[i]] <- mse(y1, r_fft_grouped$sesonal_sin)
  data$sin_err[[3]][[i]] <- mse(y2, r_fft_grouped$sesonal_cos)
  data$exp_err[[3]][[i]] <- mse(y3, r_fft_grouped$e)
  data$all_err[[3]][[i]] <- mse(y1 + y2 + y3,
                                r_fft_grouped$sesonal_sin +
                                  r_fft_grouped$sesonal_cos +
                                  r_fft_grouped$e
                                )
  
  
  data$cos_err[[4]][[i]] <- mse(y1, r_cissa_grouped$sesonal_sin)
  data$sin_err[[4]][[i]] <- mse(y2, r_cissa_grouped$sesonal_cos)
  data$exp_err[[4]][[i]] <- mse(y3, r_cissa_grouped$e)
  data$all_err[[4]][[i]] <- mse(y1 + y2 + y3,
                                r_cissa_grouped$sesonal_sin +
                                  r_cissa_grouped$sesonal_cos +
                                  r_cissa_grouped$e
                                )
  
  
  data$cos_err[[5]][[i]] <- mse(y1, r_cissa_grouped_ext$sesonal_sin)
  data$sin_err[[5]][[i]] <- mse(y2, r_cissa_grouped_ext$sesonal_cos)
  data$exp_err[[5]][[i]] <- mse(y3, r_cissa_grouped_ext$e)
  data$all_err[[5]][[i]] <- mse(y1 + y2 + y3,
                                r_cissa_grouped_ext$sesonal_sin + 
                                  r_cissa_grouped_ext$sesonal_cos +
                                  r_cissa_grouped_ext$e
                                )
  
  data$cos_err[[6]][[i]] <- mse(y1, r_fft_grouped_extended$sesonal_sin)
  data$sin_err[[6]][[i]] <- mse(y2, r_fft_grouped_extended$sesonal_cos)
  data$exp_err[[6]][[i]] <- mse(y3, r_fft_grouped_extended$e)
  data$all_err[[6]][[i]] <- mse(y1 + y2 + y3,
                                r_fft_grouped_extended$sesonal_sin +
                                  r_fft_grouped_extended$sesonal_cos +
                                  r_fft_grouped_extended$e
                                )
}
library(xtable)
# data$sin_err[[2]]
# Шаг 2: Создание примера данных
data_prev <- data
data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  exp_err = c(0, 0, 0, 0, 0, 0),
  sin_err = c(0, 0, 0, 0, 0, 0),
  cos_err = c(0, 0, 0, 0, 0, 0),
  all_err = c(0, 0, 0, 0, 0, 0)
)

data$cos_err[1] <- mean(data_prev$cos_err[[1]] |> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mean(data_prev$sin_err[[1]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$exp_err[1] <- mean(data_prev$exp_err[[1]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[1] <- mean(data_prev$all_err[[1]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[2] <- mean(data_prev$cos_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mean(data_prev$sin_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$exp_err[2] <- mean(data_prev$exp_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[2] <- mean(data_prev$all_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mean(data_prev$cos_err[[3]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mean(data_prev$sin_err[[3]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$exp_err[3] <- mean(data_prev$exp_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[3] <- mean(data_prev$all_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mean(data_prev$cos_err[[4]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mean(data_prev$sin_err[[4]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$exp_err[4] <- mean(data_prev$exp_err[[4]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[4] <- mean(data_prev$all_err[[4]]|> unlist()) |>
  formatC(format = "e", digits = 1)

data$cos_err[5] <- mean(data_prev$cos_err[[5]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mean(data_prev$sin_err[[5]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$exp_err[5] <- mean(data_prev$exp_err[[5]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[5] <- mean(data_prev$all_err[[5]]|> unlist()) |>
  formatC(format = "e", digits = 1)

data$cos_err[6] <- mean(data_prev$cos_err[[6]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mean(data_prev$sin_err[[6]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$exp_err[6] <- mean(data_prev$exp_err[[6]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[6] <- mean(data_prev$all_err[[6]]|> unlist()) |>
  formatC(format = "e", digits = 1)

for (i in 1:6){
  for (j in (i):6){
    x <- data_prev$cos_err[[i]] |> unlist()
    y <- data_prev$cos_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    if (!is.na(t_test_result$p.value) & t_test_result$p.value > 0.05){
      print(paste("cos, ", data$Метод[i], " ",
                  data$Метод[j], ", p-val = ", t_test_result$p.value))
      }
  }
}

for (i in 1:6){
  for (j in (i):6){
    x <- data_prev$sin_err[[i]] |> unlist()
    y <- data_prev$sin_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    if (!is.na(t_test_result$p.value) & t_test_result$p.value > 0.05){
      print(paste("sin, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))
      }
  }
}

[1] "sin,  SSA   SSA EOSSA , p-val =  0.171160848885882"
[1] "sin,  Fourier   CiSSA extended , p-val =  0.82766511101901"

for (i in 1:6){
  for (j in (i):6){
    x <- data_prev$exp_err[[i]] |> unlist()
    y <- data_prev$exp_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    if (!is.na(t_test_result$p.value) & t_test_result$p.value > 0.05){
      print(paste("exp, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))}
  }
}

for (i in 1:6){
  for (j in (i):6){
    x <- data_prev$all_err[[i]] |> unlist()
    y <- data_prev$all_err[[j]] |> unlist()
    t_test_result <- t.test(x, y, paired = TRUE)
    if (!is.na(t_test_result$p.value) & t_test_result$p.value > 0.05){
      print(paste("all_err, ", data$Метод[i], " ", data$Метод[j], ", p-val = ", t_test_result$p.value))
      }
  }
}


table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Fri Feb 28 17:02:44 2025
\begin{table}[ht]
\centering
\begin{tabular}{lllll}
  \hline
Метод & exp\_err & sin\_err & cos\_err & all\_err \\ 
  \hline
SSA & 5.2e-03 & 2.9e-04 & 3.6e-04 & 5.2e-03 \\ 
  SSA EOSSA & 9.5e-04 & 2.9e-04 & 3.1e-04 & 1.5e-03 \\ 
  Fourier & 1.2e-01 & 6.9e-04 & 7.2e-03 & 1.1e-01 \\ 
  CiSSA & 5.5e-02 & 1.7e-04 & 7.0e-04 & 4.6e-02 \\ 
  CiSSA extended & 2.7e-03 & 6.8e-04 & 2.1e-03 & 3.1e-03 \\ 
  Fourier extended & 3.0e-03 & 1.9e-03 & 9.6e-03 & 1.2e-02 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Пример, который не подходит ни к какому методу

n <- 96*2 + 7
L <- 89
x <- 0:(n-1)
y1 <- sin(2*pi/13 * x)
y2 <- cos(2*pi/8 * x)
y3 <- -x*x/1000
y4 <- exp(x/55)
y <- y1 + y2 + y3 + y4 
plot(x, y)

X <- hankel(y, L = L)
eps <- 1/(n+1)

s_ssa <- ssa(y, L)
# r_ssa <- reconstruct(s_ssa, groups=list(
#   sesonal_sin = c(1, 2),
#   sesonal_cos = c(3, 4)
# ))
# plot(x, r_ssa$sesonal_sin)
# plot(x, r_ssa$sesonal_cos)
# plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
e_ssa <- eossa_new(s_ssa, nested.groups = list(1:10), clust_type = "distance")
g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                       freq.bins = list(
                         sesonal_sin = c(1/12-eps, 1/12+eps),
                         sesonal_cos = c(1/8-eps, 1/8+eps)
                                        ),
                       threshold = 0.5)
r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
# plot(x, r_ssa_e$sesonal_sin)
# plot(x, r_ssa_e$sesonal_cos)



r_fft <- reconstruct_fft(x, y)
r_fft_grouped <- grouping_cissa(r_fft,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/8-eps, 1/8+eps)
                    )
                    )$t_series

# r_fft_grouped$sesonal_sin |> length()

# plot(x, r_fft_grouped$sesonal_sin)
# plot(x, r_fft_grouped$sesonal_cos)
# plot(x, r_fft_grouped$residuals)


r_cissa <- circulant_SSA(y, L)
r_cissa_grouped <- grouping_cissa(r_cissa,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/8-eps, 1/8+eps)
                    )
                    )$t_series


r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/8-eps, 1/8+eps)
                    )
                    )$t_series

# plot(x, r_cissa_grouped$sesonal_sin)
# plot(x, r_cissa_grouped$sesonal_cos)
# plot(x, r_cissa_grouped$residuals)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA EOSSA","Fourier", "CiSSA", "CiSSA extended"),
  sin_err = c (20, 20, 1, 1),
  cos_err = c(1, 1, 1, 1)
)

data$cos_err[1] <- mse(y1, r_ssa_e$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mse(y2, r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)


data$cos_err[2] <- mse(y1, r_fft_grouped$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(y2, r_fft_grouped$sesonal_cos) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mse(y1, r_cissa_grouped$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mse(y2, r_cissa_grouped$sesonal_cos) |> 
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mse(y1, r_cissa_grouped_ext$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mse(y2, r_cissa_grouped_ext$sesonal_cos) |> 
  formatC(format = "e", digits = 1)



table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Fri Feb 28 17:02:45 2025
\begin{table}[ht]
\centering
\begin{tabular}{lll}
  \hline
Метод & sin\_err & cos\_err \\ 
  \hline
SSA EOSSA & 2.1e-21 & 7.4e-22 \\ 
  Fourier & 1.6e-02 & 4.2e-01 \\ 
  CiSSA & 2.2e-02 & 3.7e-02 \\ 
  CiSSA extended & 2.9e-03 & 4.5e-03 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Точная разделимость. Примеры

library(Rssa)

do_one_experiment <- function(methods_list, methods_names){
  
}

n <- 96*2
L <- 96
x <- 0:(n-1)
y1 <- sin(2*pi/12 * x)
y2 <- cos(2*pi/3 * x)/2
y <- y1 + y2
X <- hankel(y, L = L)
eps <- 1/(n+1)

s_ssa <- ssa(y[1:(n-1)], L)
r_ssa <- reconstruct(s_ssa, groups = list(
  F1 = c(1, 2),
  F2 = c(3, 4)
))
# r_ssa <- reconstruct(s_ssa, groups=list(
#   sesonal_sin = c(1, 2),
#   sesonal_cos = c(3, 4)
# ))
# plot(x, r_ssa$F1, type="l")
# plot(x, r_ssa$F2, type= "l")
# plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
e_ssa <- eossa_new(s_ssa, nested.groups = list(1:4), clust_type = "distance")
g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                       freq.bins = list(
                         sesonal_sin = c(1/12-eps, 1/12+eps),
                         sesonal_cos = c(1/3-eps, 1/3+eps)
                                        ),
                       threshold = 0.5)
r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
# plot(x, r_ssa_e$sesonal_sin)
# plot(x, r_ssa_e$sesonal_cos)



r_fft <- reconstruct_fft(x, y)
r_fft_grouped <- grouping_cissa(r_fft,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

r_fft_extended <- reconstruct_fft(x, y, TRUE)
r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# r_fft_grouped$sesonal_sin |> length()

# plot(x, r_fft_grouped$sesonal_sin)
# plot(x, r_fft_grouped$sesonal_cos)
# plot(x, r_fft_grouped$residuals)


r_cissa <- circulant_SSA(y, L)
r_cissa_grouped <- grouping_cissa(r_cissa,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series


r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series






# plot(x, r_cissa_grouped$sesonal_sin)
# plot(x, r_cissa_grouped$sesonal_cos)
# plot(x, r_cissa_grouped$residuals)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA, \nLw, Kw in N", "SSA EOSSA, \nLw, Kw in N","Fourier, Nw in N", "CiSSA, Lw in N", "CiSSA extended, Lw in N", "Fourier extended, Nw in N",
            "SSA, \nLw in N, Kw not in N", "SSA EOSSA, \nLw in N, Kw not in N","Fourier, Nw not in N", "CiSSA, Lw not in N", "CiSSA extended, Lw not in N", "Fourier extended, Nw not in N"),
  sin_err = c (1, 20, 20, 1, 1, 1, 1, 20, 20, 1, 1, 1),
  cos_err = c(1, 1, 1, 1, 1, 1, 1, 20, 20, 1, 1, 1),
  all_err = c(1, 1, 1, 1, 1, 1, 1, 20, 20, 1, 1, 1)
)

data$sin_err[1] <- mse(y1[1:(n-1)], r_ssa$F1) |>
  formatC(format = "e", digits = 1) 
data$cos_err[1] <- mse(y2[1:(n-1)], r_ssa$F2) |>
  formatC(format = "e", digits = 1)
data$all_err[1] <- mse(y1[1:(n-1)] + y2[1:(n-1)], r_ssa$F1 + r_ssa$F2) |>
  formatC(format = "e", digits = 1)

data$cos_err[2] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$all_err[2] <- mse(y1[1:(n-1)] + y2[1:(n-1)],
                       r_ssa_e$sesonal_sin + r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mse(y1, r_fft_grouped$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mse(y2, r_fft_grouped$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$all_err[3] <- mse(y1 + y2,
                       r_fft_grouped$sesonal_sin + r_fft_grouped$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mse(y1, r_cissa_grouped$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mse(y2, r_cissa_grouped$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$all_err[4] <- mse(y1 + y2,
                       r_cissa_grouped$sesonal_sin + r_cissa_grouped$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[5] <- mse(y1, r_cissa_grouped_ext$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mse(y2, r_cissa_grouped_ext$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$all_err[5] <- mse(y1 + y2,
                       r_cissa_grouped_ext$sesonal_sin + r_cissa_grouped_ext$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[6] <- mse(y1, r_fft_grouped_extended$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mse(y2, r_fft_grouped_extended$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$all_err[6] <- mse(y1 + y2,
                       r_fft_grouped_extended$sesonal_sin + r_fft_grouped_extended$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)



#------------------

eps <- eps*2

s_ssa <- ssa(y, L)
r_ssa <- reconstruct(s_ssa, groups = list(
  F1 = c(1, 2),
  F2 = c(3, 4)
))
# r_ssa <- reconstruct(s_ssa, groups=list(
#   sesonal_sin = c(1, 2),
#   sesonal_cos = c(3, 4)
# ))
# plot(x, r_ssa$F1, type="l")
# plot(x, r_ssa$F2, type= "l")
# plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
e_ssa <- eossa_new(s_ssa, nested.groups = list(1:4), clust_type = "distance")
g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                       freq.bins = list(
                         sesonal_sin = c(1/12-eps, 1/12+eps),
                         sesonal_cos = c(1/3-eps, 1/3+eps)
                                        ),
                       threshold = 0.5)
r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
# plot(x, r_ssa_e$sesonal_sin)
# plot(x, r_ssa_e$sesonal_cos)



r_fft <- reconstruct_fft(x[1:(n-1)], y[1:(n-1)])
r_fft_grouped <- grouping_cissa(r_fft,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

r_fft_extended <- reconstruct_fft(x[1:(n-1)], y[1:(n-1)], TRUE)
r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# r_fft_grouped$sesonal_sin |> length()

# plot(x, r_fft_grouped$sesonal_sin)
# plot(x, r_fft_grouped$sesonal_cos)
# plot(x, r_fft_grouped$residuals)


r_cissa <- circulant_SSA(y, L+1)
r_cissa_grouped <- grouping_cissa(r_cissa,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series


r_cissa_ext <- circulant_SSA(y, L+1, extend_flag = TRUE)
r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                    groups = list(
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series





data$sin_err[7] <- mse(y1, r_ssa$F1) |>
  formatC(format = "e", digits = 1) 
data$cos_err[7] <- mse(y2, r_ssa$F2) |>
  formatC(format = "e", digits = 1)
data$all_err[7] <- mse(y1 + y2, r_ssa$F1 + r_ssa$F2) |>
  formatC(format = "e", digits = 1)

data$cos_err[8] <- mse(y1, r_ssa_e$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[8] <- mse(y2, r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$all_err[8] <- mse(y1 + y2,
                       r_ssa_e$sesonal_sin + r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)


data$cos_err[9] <- mse(y1[1:(n-1)], r_fft_grouped$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[9] <- mse(y2[1:(n-1)], r_fft_grouped$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$all_err[9] <- mse(y1[1:(n-1)] + y2[1:(n-1)],
                       r_fft_grouped$sesonal_sin + r_fft_grouped$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[10] <- mse(y1, r_cissa_grouped$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[10] <- mse(y2, r_cissa_grouped$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$all_err[10] <- mse(y1 + y2,
                       r_cissa_grouped$sesonal_sin + r_cissa_grouped$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[11] <- mse(y1, r_cissa_grouped_ext$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[11] <- mse(y2, r_cissa_grouped_ext$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$all_err[11] <- mse(y1 + y2,
                       r_cissa_grouped_ext$sesonal_sin + r_cissa_grouped_ext$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[12] <- mse(y1[1:(n-1)], r_fft_grouped_extended$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[12] <- mse(y2[1:(n-1)], r_fft_grouped_extended$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$all_err[12] <- mse(y1[1:(n-1)] + y2[1:(n-1)],
                       r_fft_grouped_extended$sesonal_sin + r_fft_grouped_extended$sesonal_cos
                       ) |>
  formatC(format = "e", digits = 1)






table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Fri Feb 28 17:02:48 2025
\begin{table}[ht]
\centering
\begin{tabular}{llll}
  \hline
Метод & sin\_err & cos\_err & all\_err \\ 
  \hline
SSA, 
Lw, Kw in N & 6.8e-30 & 1.5e-29 & 1.8e-29 \\ 
  SSA EOSSA, 
Lw, Kw in N & 1.5e-29 & 7.5e-30 & 2.0e-29 \\ 
  Fourier, Nw in N & 1.7e-28 & 3.5e-28 & 5.1e-28 \\ 
  CiSSA, Lw in N & 1.9e-29 & 5.3e-30 & 2.1e-29 \\ 
  CiSSA extended, Lw in N & 2.0e-04 & 8.6e-04 & 1.1e-03 \\ 
  Fourier extended, Nw in N & 6.2e-04 & 2.6e-03 & 3.2e-03 \\ 
  SSA, 
Lw in N, Kw not in N & 2.2e-06 & 2.2e-06 & 2.0e-29 \\ 
  SSA EOSSA, 
Lw in N, Kw not in N & 1.5e-29 & 8.8e-30 & 1.9e-29 \\ 
  Fourier, Nw not in N & 9.4e-03 & 3.5e-03 & 1.3e-02 \\ 
  CiSSA, Lw not in N & 1.7e-02 & 7.0e-03 & 2.3e-02 \\ 
  CiSSA extended, Lw not in N & 2.4e-03 & 6.9e-04 & 3.1e-03 \\ 
  Fourier extended, Nw not in N & 1.1e-05 & 4.9e-04 & 4.9e-04 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Те же примеры с шумом

data <- list(
  Метод = c("SSA","SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = list(list(), list(), list(), list(), list(), list()),
  cos_err = list(list(), list(), list(), list(), list(), list()),
  all_err = list(list(), list(), list(), list(), list(), list())
)

for (i in 1:100){
  set.seed(i)
  
  n <- 96*2
  L <- 96
  x <- 0:(n-1)
  y1 <- sin(2*pi/12 * x)
  y2 <- cos(2*pi/3 * x)/2
  y <- y1 + y2 + rnorm(n, 0, 0.1)
  X <- hankel(y, L = L)
  eps <- 1/(n+1)
  
  s_ssa <- ssa(y[1:(n-1)], L)
  r_ssa <- reconstruct(s_ssa, groups = list(
    F1 = c(1, 2),
    F2 = c(3, 4)
  ))
  # r_ssa <- reconstruct(s_ssa, groups=list(
  #   sesonal_sin = c(1, 2),
  #   sesonal_cos = c(3, 4)
  # ))
  # plot(x[1:(n-1)], r_ssa$F1, type = "l")
  # plot(x[1:(n-1)], r_ssa$F2, type = "l")
  # plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
  e_ssa <- eossa_new(s_ssa, nested.groups = list(1:4), clust_type = "distance")
  g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                         freq.bins = list(
                           sesonal_sin = c(1/12-eps, 1/12+eps),
                           sesonal_cos = c(1/3-eps, 1/3+eps)
                                          ),
                         threshold = 0.5)
  r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
  # plot(x, r_ssa_e$sesonal_sin)
  # plot(x, r_ssa_e$sesonal_cos)
  
  
  
  r_fft <- reconstruct_fft(x, y)
  r_fft_grouped <- grouping_cissa(r_fft,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  r_fft_extended <- reconstruct_fft(x, y, TRUE)
  r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  # r_fft_grouped$sesonal_sin |> length()
  
  # plot(x, r_fft_grouped$sesonal_sin)
  # plot(x, r_fft_grouped$sesonal_cos)
  # plot(x, r_fft_grouped$residuals)
  
  
  r_cissa <- circulant_SSA(y, L)
  r_cissa_grouped <- grouping_cissa(r_cissa,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  
  
  r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
  r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                      groups = list(
                        sesonal_sin = c(1/12-eps, 1/12+eps),
                        sesonal_cos = c(1/3-eps, 1/3+eps)
                      )
                      )$t_series
  data$cos_err[[1]][[i]] <- 
    min(mse(y1[1:(n-1)], r_ssa$F1), mse(y1[1:(n-1)], r_ssa$F2))
  data$sin_err[[1]][[i]] <-
    min(mse(y2[1:(n-1)], r_ssa$F1), mse(y2[1:(n-1)], r_ssa$F2))
  data$all_err[[1]][[i]] <-
    min(mse(y1[1:(n-1)]+y2[1:(n-1)], r_ssa$F1+ r_ssa$F2))
  
  data$cos_err[[2]][[i]] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin)
  data$sin_err[[2]][[i]] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos)
  data$all_err[[2]][[i]] <- mse(y1[1:(n-1)]+ y2[1:(n-1)],
                                r_ssa_e$sesonal_sin + r_ssa_e$sesonal_cos)
  
  
  data$cos_err[[3]][[i]] <- mse(y1, r_fft_grouped$sesonal_sin)
  data$sin_err[[3]][[i]] <- mse(y2, r_fft_grouped$sesonal_cos)
  data$all_err[[3]][[i]] <- mse(y1 + y2,
                                r_fft_grouped$sesonal_cos +  r_fft_grouped$sesonal_sin)
  
  
  data$cos_err[[4]][[i]] <- mse(y1, r_cissa_grouped$sesonal_sin)
  data$sin_err[[4]][[i]] <- mse(y2, r_cissa_grouped$sesonal_cos)
  data$all_err[[4]][[i]] <- mse(y1 + y2,
                          r_cissa_grouped$sesonal_sin +  r_cissa_grouped$sesonal_cos
                                )
  
  data$cos_err[[5]][[i]] <- mse(y1, r_cissa_grouped_ext$sesonal_sin)
  data$sin_err[[5]][[i]] <- mse(y2, r_cissa_grouped_ext$sesonal_cos)
  data$all_err[[5]][[i]] <- mse(y1 + y2,
                          r_cissa_grouped_ext$sesonal_sin + r_cissa_grouped_ext$sesonal_cos
                                )
  
  data$cos_err[[6]][[i]] <- mse(y1, r_fft_grouped_extended$sesonal_sin)
  data$sin_err[[6]][[i]] <- mse(y2, r_fft_grouped_extended$sesonal_cos)
  data$all_err[[6]][[i]] <- mse(y1 + y2,
                          r_fft_grouped_extended$sesonal_sin + r_fft_grouped_extended$sesonal_cos
                                )
}
library(xtable)

# Шаг 2: Создание примера данных
data_prev <- data

data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  sin_err = c(0, 0, 0, 0, 0, 0),
  cos_err = c(0, 0, 0, 0, 0, 0),
  all_err = c(0, 0, 0, 0, 0, 0)
)

data$cos_err[1] <- mean(data_prev$cos_err[[1]] |> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mean(data_prev$sin_err[[1]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[1] <- mean(data_prev$all_err[[1]]|> unlist()) |>
  formatC(format = "e", digits = 1)

data$cos_err[2] <- mean(data_prev$cos_err[[2]] |> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mean(data_prev$sin_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[2] <- mean(data_prev$all_err[[2]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mean(data_prev$cos_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mean(data_prev$sin_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)
data$all_err[3] <- mean(data_prev$all_err[[3]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mean(data_prev$cos_err[[4]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mean(data_prev$sin_err[[4]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$all_err[4] <- mean(data_prev$all_err[[4]]|> unlist()) |>
  formatC(format = "e", digits = 1)


data$cos_err[5] <- mean(data_prev$cos_err[[5]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mean(data_prev$sin_err[[5]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$all_err[5] <- mean(data_prev$all_err[[5]]|> unlist()) |>
  formatC(format = "e", digits = 1)

data$cos_err[6] <- mean(data_prev$cos_err[[6]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mean(data_prev$sin_err[[6]]|> unlist()) |> 
  formatC(format = "e", digits = 1)
data$all_err[6] <- mean(data_prev$all_err[[6]]|> unlist()) |>
  formatC(format = "e", digits = 1)

Ассимптотическая разделимость (бесполезно)

n <- 96*2
x <- 0:(n-1)
L <- 96
y1 <- sin(2*pi/12 * x)
y2 <- cos(2*pi/3 * x)/2
y3 <- exp(x/10) + 1 
y <- y1 + y2 + y3
eps <- 1/(n+1)

s_ssa <- ssa(y[1:(n-1)], L)
r_ssa <- reconstruct(s_ssa, groups=list(
  e = c(1),
  sesonal_sin = c(2, 3),
  sesonal_cos = c(4, 5)
))
# plot(x[1:(n-1)], r_ssa$e, type = "l")
# plot(x[1:(n-1)], r_ssa$sesonal_sin, type = "l")
# plot(x[1:(n-1)], r_ssa$sesonal_cos, type = "l")
# plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
e_ssa <- eossa_new(s_ssa, nested.groups = list(1:7), clust_type = "distance")
g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                       freq.bins = list(
                         e = c(0, 1/12-eps-eps),
                         sesonal_sin = c(1/12-eps, 1/12+eps),
                         sesonal_cos = c(1/3-eps, 1/3+eps)
                                        ),
                       threshold = 0.5)
r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
# plot(x, r_ssa_e$sesonal_sin)
# plot(x, r_ssa_e$sesonal_cos)



r_fft <- reconstruct_fft(x, y)
r_fft_grouped <- grouping_cissa(r_fft,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

r_fft_extended <- reconstruct_fft(x, y, TRUE)
r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# r_fft_grouped$sesonal_sin |> length()

# plot(x, r_fft_grouped$sesonal_sin)
# plot(x, r_fft_grouped$sesonal_cos)
# plot(x, r_fft_grouped$residuals)


r_cissa <- circulant_SSA(y, L)
r_cissa_grouped <- grouping_cissa(r_cissa,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series


r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# plot(x, r_cissa_grouped$sesonal_sin)
# plot(x, r_cissa_grouped$sesonal_cos)
# plot(x, r_cissa_grouped$residuals)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  e_err = c(1, 1, 1, 1, 1, 1),
  sin_err = c (1, 20, 20, 1, 1, 1),
  cos_err = c(1, 1, 1, 1, 1, 1),
  all_err = c(1, 1, 1, 1, 1, 1)
)

data$cos_err[1] <- mse(y1[1:(n-1)], r_ssa$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mse(y2[1:(n-1)], r_ssa$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[1] <- mse(y3[1:(n-1)], r_ssa$e) |>
  formatC(format = "e", digits = 1)
data$all_err[1] <- mse(y1[1:(n-1)]+ y2[1:(n-1)] + y3[1:(n-1)],
                     r_ssa$sesonal_sin + r_ssa$sesonal_cos + r_ssa$e) |>
  formatC(format = "e", digits = 1)

data$cos_err[2] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[2] <- mse(y3[1:(n-1)], r_ssa_e$e) |>
  formatC(format = "e", digits = 1)
data$all_err[2] <- mse(y1[1:(n-1)]+ y2[1:(n-1)] + y3[1:(n-1)],
                     r_ssa_e$sesonal_sin + r_ssa_e$sesonal_cos + r_ssa_e$e
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mse(y1, r_fft_grouped$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mse(y2, r_fft_grouped$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[3] <- mse(y3, r_fft_grouped$e) |>
  formatC(format = "e", digits = 1)
data$all_err[3] <- mse(y1 + y2 + y3, 
                       r_fft_grouped$sesonal_sin + r_fft_grouped$sesonal_cos + r_fft_grouped$e) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mse(y1, r_cissa_grouped$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mse(y2, r_cissa_grouped$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[4] <- mse(y3, r_cissa_grouped$e) |>
  formatC(format = "e", digits = 1)
data$all_err[4] <- mse(y1 + y2 + y3, 
                       r_cissa_grouped$sesonal_sin + 
                         r_cissa_grouped$sesonal_cos + 
                         r_cissa_grouped$e
                       ) |>
  formatC(format = "e", digits = 1)

data$cos_err[5] <- mse(y1, r_cissa_grouped_ext$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mse(y2, r_cissa_grouped_ext$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[5] <- mse(y3, r_cissa_grouped_ext$e) |>
  formatC(format = "e", digits = 1)
data$all_err[5] <- mse(y1 + y2 + y3, 
                       r_cissa_grouped_ext$sesonal_sin + 
                         r_cissa_grouped_ext$sesonal_cos +
                          r_cissa_grouped_ext$e
                       ) |>
  formatC(format = "e", digits = 1)

data$cos_err[6] <- mse(y1, r_fft_grouped_extended$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mse(y2, r_fft_grouped_extended$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[6] <- mse(y3, r_fft_grouped_extended$e) |>
  formatC(format = "e", digits = 1)
data$all_err[6] <- mse(y1 + y2 + y3, 
                       r_fft_grouped_extended$sesonal_sin +
                         r_fft_grouped_extended$sesonal_cos +
                         r_fft_grouped_extended$e
                       ) |>
  formatC(format = "e", digits = 1)



table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)

% latex table generated in R 4.2.2 by xtable 1.8-4 package
% Tue Feb 25 11:51:22 2025
\begin{table}[ht]
\centering
\begin{tabular}{lllll}
  \hline
Метод & e\_err & sin\_err & cos\_err & all\_err \\ 
  \hline
SSA & 5.3e-01 & 1.5e-01 & 6.6e-01 & 3.8e-02 \\ 
  SSA EOSSA & 8.0e-05 & 5.8e-04 & 2.1e-04 & 7.3e-04 \\ 
  Fourier & 1.4e+14 & 7.8e+11 & 8.4e+12 & 1.4e+14 \\ 
  CiSSA & 1.0e+14 & 3.3e+10 & 8.5e+11 & 8.3e+13 \\ 
  CiSSA extended & 2.5e+12 & 3.0e+08 & 6.2e+10 & 1.8e+12 \\ 
  Fourier extended & 3.0e+12 & 2.9e+11 & 3.4e+12 & 6.6e+12 \\ 
   \hline
\end{tabular}
\caption{Example Table} 
\end{table}

Несколько непериодик

n <- 96*2
x <- 0:(n-1)
L <- 96
y1 <- sin(2*pi/12 * x)
y2 <- cos(2*pi/3 * x)/2
y3 <- exp(x/100) + 1
y <- y1 + y2 + y3
eps <- 1/(n+1)

s_ssa <- ssa(y[1:(n-1)], L)
r_ssa <- reconstruct(s_ssa, groups=list(
  e = c(1),
  sesonal_sin = c(2, 3),
  sesonal_cos = c(4, 5)
))
# plot(x[1:(n-1)], r_ssa$e, type = "l")
# plot(x[1:(n-1)], r_ssa$sesonal_sin, type = "l")
# plot(x[1:(n-1)], r_ssa$sesonal_cos, type = "l")
# plot(wcor(s_ssa, groups = 1:10), scales = list(at = c(10, 20, 30)))
e_ssa <- eossa_new(s_ssa, nested.groups = list(1:10), clust_type = "distance")
g_sesonal_e <- grouping.auto(e_ssa, base = "eigen",
                       freq.bins = list(
                         e = c(0, 1/12-eps-eps),
                         sesonal_sin = c(1/12-eps, 1/12+eps),
                         sesonal_cos = c(1/3-eps, 1/3+eps)
                                        ),
                       threshold = 0.5)
r_ssa_e <- reconstruct(e_ssa, groups=g_sesonal_e)
r_ssa_e
# plot(x, r_ssa_e$sesonal_sin)
# plot(x, r_ssa_e$sesonal_cos)



r_fft <- reconstruct_fft(x, y)
r_fft_grouped <- grouping_cissa(r_fft,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

r_fft_extended <- reconstruct_fft(x, y, TRUE)
r_fft_grouped_extended <- grouping_cissa(r_fft_extended,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# r_fft_grouped$sesonal_sin |> length()

# plot(x, r_fft_grouped$sesonal_sin)
# plot(x, r_fft_grouped$sesonal_cos)
# plot(x, r_fft_grouped$residuals)


r_cissa <- circulant_SSA(y, L)
r_cissa_grouped <- grouping_cissa(r_cissa,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series


r_cissa_ext <- circulant_SSA(y, L, extend_flag = TRUE)
r_cissa_grouped_ext <- grouping_cissa(r_cissa_ext,
                    groups = list(
                      e = c(0, 1/12-eps-eps),
                      sesonal_sin = c(1/12-eps, 1/12+eps),
                      sesonal_cos = c(1/3-eps, 1/3+eps)
                    )
                    )$t_series

# plot(x, r_cissa_grouped$sesonal_sin)
# plot(x, r_cissa_grouped$sesonal_cos)
# plot(x, r_cissa_grouped$residuals)

library(xtable)

# Шаг 2: Создание примера данных
data <- data.frame(
  Метод = c("SSA", "SSA EOSSA","Fourier", "CiSSA", "CiSSA extended", "Fourier extended"),
  e_err = c(1, 1, 1, 1, 1, 1),
  sin_err = c (1, 20, 20, 1, 1, 1),
  cos_err = c(1, 1, 1, 1, 1, 1),
  all_err = c(1, 1, 1, 1, 1, 1)
)

data$cos_err[1] <- mse(y1[1:(n-1)], r_ssa$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[1] <- mse(y2[1:(n-1)], r_ssa$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[1] <- mse(y3[1:(n-1)], r_ssa$e) |>
  formatC(format = "e", digits = 1)
data$all_err[1] <- mse(y1[1:(n-1)]+ y2[1:(n-1)] + y3[1:(n-1)],
                     r_ssa$sesonal_sin + r_ssa$sesonal_cos + r_ssa$e) |>
  formatC(format = "e", digits = 1)

data$cos_err[2] <- mse(y1[1:(n-1)], r_ssa_e$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[2] <- mse(y2[1:(n-1)], r_ssa_e$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[2] <- mse(y3[1:(n-1)], r_ssa_e$e) |>
  formatC(format = "e", digits = 1)
data$all_err[2] <- mse(y1[1:(n-1)]+ y2[1:(n-1)] + y3[1:(n-1)],
                     r_ssa_e$sesonal_sin + r_ssa_e$sesonal_cos + r_ssa_e$e
                       ) |>
  formatC(format = "e", digits = 1)


data$cos_err[3] <- mse(y1, r_fft_grouped$sesonal_sin) |>
  formatC(format = "e", digits = 1)
data$sin_err[3] <- mse(y2, r_fft_grouped$sesonal_cos) |>
  formatC(format = "e", digits = 1)
data$e_err[3] <- mse(y3, r_fft_grouped$e) |>
  formatC(format = "e", digits = 1)
data$all_err[3] <- mse(y1 + y2 + y3, 
                       r_fft_grouped$sesonal_sin + r_fft_grouped$sesonal_cos + r_fft_grouped$e) |>
  formatC(format = "e", digits = 1)


data$cos_err[4] <- mse(y1, r_cissa_grouped$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[4] <- mse(y2, r_cissa_grouped$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[4] <- mse(y3, r_cissa_grouped$e) |>
  formatC(format = "e", digits = 1)
data$all_err[4] <- mse(y1 + y2 + y3, 
                       r_cissa_grouped$sesonal_sin + 
                         r_cissa_grouped$sesonal_cos + 
                         r_cissa_grouped$e
                       ) |>
  formatC(format = "e", digits = 1)

data$cos_err[5] <- mse(y1, r_cissa_grouped_ext$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[5] <- mse(y2, r_cissa_grouped_ext$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[5] <- mse(y3, r_cissa_grouped_ext$e) |>
  formatC(format = "e", digits = 1)
data$all_err[5] <- mse(y1 + y2 + y3, 
                       r_cissa_grouped_ext$sesonal_sin + 
                         r_cissa_grouped_ext$sesonal_cos +
                          r_cissa_grouped_ext$e
                       ) |>
  formatC(format = "e", digits = 1)

data$cos_err[6] <- mse(y1, r_fft_grouped_extended$sesonal_sin) |> 
  formatC(format = "e", digits = 1)
data$sin_err[6] <- mse(y2, r_fft_grouped_extended$sesonal_cos) |> 
  formatC(format = "e", digits = 1)
data$e_err[6] <- mse(y3, r_fft_grouped_extended$e) |>
  formatC(format = "e", digits = 1)
data$all_err[6] <- mse(y1 + y2 + y3, 
                       r_fft_grouped_extended$sesonal_sin +
                         r_fft_grouped_extended$sesonal_cos +
                         r_fft_grouped_extended$e
                       ) |>
  formatC(format = "e", digits = 1)



table_latex <- xtable(data, caption = "Example Table")

# Шаг 4: Вывод таблицы в LaTeX файл
print(table_latex, include.rownames = FALSE)